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:41 MST