Re: Downsizing a ncdf file

From: Dennis Shea (shea AT cgd.ucar.edu)
Date: Mon Feb 14 2005 - 14:21:09 MST

  • Next message: Saji Hameed: "Color for _FillValue"

    Hi Adolf

    In addition to the response by Saji N. Hameed, Matt Stumbaugh
    sent to me the attached script and pdf file.

    THX to both Saji and Matt!

    Regards
    Dennis Shea

    Dennis,

    I'm not sure if this could help Adolf but reading his question it seemed
    part of his problem was taking higher resolution data and converting it
    to lower (15x15 --> 5x5). Basically, it takes arguments of polygon data
    (latitude, longitude) of unknown resolution, coordinate arrays (lat,
    lon) for a known grid, and sorts the polygons into a grid corresponding
    to the coordinate arrays(lat x lon). It is just a way of taking fairly
    large polygon data sets and reducing them to a coarser grid so I can see
    where data is available without having to dig into the higher resolution
    data. If you think it will help please pass this on. I've attached a
    plot. In the figure, it appears like I've just used polyline and
    polymarkers on the same data but in fact I used them respectively on the
    polygon data and the polygon data devolved raster data. I have to get
    other work done now but if this code is of interest I can clean it up.

    Sorry to hear about Sylvia. She will be missed but I'm sure she'll do
    well wherever she lands.

    Matt

    ;---------------------------------------------------------------;
    ; poly2layer ;
    ;---------------------------------------------------------------;
    ; When positioned in a time loop, "poly2layer()" populates
    ; a cube (t X lon X lat) of known resolution by layer (lon X lat)
    ; with data from longitude & latitude arrays of unknown resolution
    ; with rows of vertice position data & columns of consecutive
    ; polygons indexed daily in a vector.
    ; There is no interpolation between data that falls in a
    ; discretized subcube of resolution (dlon, dlat).
    ;
    ;undef("poly2layer")
    function poly2layer(yPoly[*][*]:float, xPoly[*][*]:float, \
        latArray[*]:float, lonArray[*]:float, polyTime[*]:integer, \
        polyNumber[*]:integer, cubeTime:integer)
       
    local iiPolyTime, polyXx, polyYy, nLayerRow, latBottom, latTop, \
        layerRow, nPolyTime, cPolyYy, cPolyXx, icPolyYy, scPolyXx, \
        nLayerCol, lonE, lonW

    begin
        layer = new((/dimsizes(lonArray),dimsizes(latArray)/),"float")
        iiPolyTime = ind(polyTime.eq.cubeTime)
        if any(ismissing(iiPolyTime).eq.False).eq.True
            ;print("Been here")
            polyXx = xPoly(:,iiPolyTime)
            polyYy = yPoly(:,iiPolyTime)
            polyNa = polyNumber(iiPolyTime) ; polyNumbers from day
           
            do nLayerRow = 0, dimsizes(latArray)-2
                latBottom= latArray(nLayerRow)
                latTop = latArray(nLayerRow+1)
                ; This is to center vertices around the devolved grid
                ; vertices...
                latD = (latTop - latBottom)*0.5
                latBottom= latBottom - latD
                latTop = latTop - latD
                layerRow = new((/dimsizes(lonArray)/),"integer")
           
                do nPolyTime = 0, dimsizes(iiPolyTime)-1
                    cPolyYy = polyYy(:,nPolyTime)
                    cPolyXx = polyXx(:,nPolyTime)
                    icPolyYy =
    ind(cPolyYy.lt.latTop.and.cPolyYy.ge.latBottom)
                    polyN = polyNa(nPolyTime) ; polyNumber within day
                   
                    if any(ismissing(icPolyYy).eq.False).eq.True
                        scPolyXx = cPolyXx(icPolyYy)
                        do nLayerCol = 0, dimsizes(lonArray)-2
                            lonE = lonArray(nLayerCol+1)
                            lonW = lonArray(nLayerCol)
                            ; This is to center vertices around the devolved
    grid
                            ; vertices...
                            lonD = (lonE-lonW)*0.5
                            lonW = lonW - lonD
                            lonE = lonE - lonD
                            do nscPolyXx = 0, dimsizes(scPolyXx)-1
                                sscPolyXx = scPolyXx(nscPolyXx)
                                if((sscPolyXx.ge.lonW).and.(sscPolyXx.lt.lonE))
                                    layerRow(nLayerCol) = polyN ; polyNumber
    as data
                                    print("Bingo!")
                                    ; A function to fill between values is
    needed.
                                end if
                            end do
                        end do
                        delete(scPolyXx)
                    end if
                    delete(cPolyYy)
                    delete(cPolyXx)
                    delete(icPolyYy)
                end do
                layer(:,nLayerRow) = layerRow
                delete(layerRow)
            end do
        end if
        return(layer)
    end

    Dennis Shea wrote:

    >>I have model output ncdf files (lat,lon,time,levels and parameters),
    >>on say 5x5 minute grid (limited domain).
    >>I want to calculate area mean values for eg. 15x15 grid.
    >>I know ncks, extracts for example every third lat x lon,
    >>but I want the mean values. So I didn't find a way
    >>to do it with nco.
    >>(E.g extracting 3 times with ncks (shifted start) and
    >>averaging with ncea the 3 files, gives the mean values, but
    >>the coordinates are not averaged, further the handling of missing values
    >>creates another problem)
    >>
    >>Does anybody know code from NCL or..
    >>which would perform just this and may be general enough to
    >>be applied on rather arbitrary (my) ncdf files.
    >>
    >
    >Umm, I think that it is possible to use the NCO to accomplish
    >this task. I will ask someone who does this type of thing.
    >
    >---
    >
    >As far as NCL is concerned, there are numerous
    >approaches depending upon what your objective is.
    >
    >Let x(time,lev,lat,lon) where each dimension is a classical
    >netCDF coordinate variable [ie: one-dimensional and monotonically
    >{in/de}creasing]. If missing values are present the
    >_FillValue attribute should be associated with the variable x.
    >
    >
    >see: http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/wareaave.html
    >
    >
    >load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    >
    > latS = 20. ; 15 x 15 area
    > latN = 35.
    > lonL = 190.
    > lonR = 205.
    > opt = ; 0 or 1 (see write-up of wgt_areaave)
    >
    > f = addfile ("foo.nc" , "r")
    > x = f->SomeVariable(:,:,{latS:latN},{lonL:lonR} ; (time,lev,lat,lon)
    >
    > lat = x&lat ; extract latitudes associated with x
    >
    > clat = cos(lat*0.01745329) ; cos latitude wgting
    >
    > xAreaAve = wgt_areaave_Wrap(x, clat, 1., opt) ; ==> (time,lev)
    > printVarSummary(xAreaAve)
    >
    >----
    >If your lat/lon are 2D see
    >
    >http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/wareaave2.html
    >
    >----
    >Cheers!
    >
    >_______________________________________________
    >ncl-talk mailing list
    >ncl-talk@ucar.edu
    >http://mailman.ucar.edu/mailman/listinfo/ncl-talk
    >




    _______________________________________________
    ncl-talk mailing list
    ncl-talk@ucar.edu
    http://mailman.ucar.edu/mailman/listinfo/ncl-talk



    This archive was generated by hypermail 2b29 : Mon Feb 14 2005 - 15:29:44 MST