*** region_ind works just fine ***
I have attached a sample script.
Entire grid
(0) latitude: min=20.7095 max=69.9062
(0) longitude: min=203.718 max=322.936
Subset region
(0) latitude: min=36.1247 max=47.7835
(0) longitude: min=268.495 max=281.103
Sample output:
http://www.cgd.ucar.edu/~shea/MM5-0.png
http://www.cgd.ucar.edu/~shea/MM5-1.png
[...] wrote:
> Sho,
>
> There indeed seems to be a problem with the region_ind function. The
> current version requires the subscripts on lat and lon to be in (yc, xc)
> order, if I am not mistaken. Your arrays are in (xc, yc) order. This
> constraint seems to be missing from the current documentation.
>
> It turns out that the function works symmetrically. So I think you will
> get correct results if you transpose the function output array, and
> assign i1/i2 and j1/j2 differently, like this:
>
> i1 = ji(0) ; lon start
> i2 = ji(1) ; lon last
>
> j1 = ji(2) ; lat start
> j2 = ji(3) ; lat last
>
> Please consider this a temporary fix until the NCL team decides what to
> do with this function. They might decide to just document the sequence
> association between input and output subscripts, which would make the
> above code legal.
>
> Also, please try to write to me only on the users list, unless a private
> message is necessary. I do not work for Unidata, I am just a regular
> NCL user like yourself. Thank you.
>
> --[...]
>
> Sho Kawazoe wrote:
>> [...],
>>
>> Thank you for your continued help on this issue. It looks like I'm
>> now on
>> the right track, but I do have one more question that perhaps you
>> could help
>> me out on.
>>
>> Your example (as well as the example online) seem to be useful if the
>> NetCDF
>> file has grid points of standard lat/lon coordinate (that is, ranging
>> from
>> -90 to 90 in the lat, and -180 to 180 for lon). When I apply the same
>> method, the grid point I am cutting out (using specified lat points of
>> 37-47
>> and lon points of 269-279), seems to read a 10degree by 10 degree box
>> somewhere around the Montana/Alberta area. The min and max of the
>> converted
>> lat/lon variable is attached, but it looks like the script is doing some
>> wierd lat/lon conversion. I've attached the script as well as the output.
>> I've fiddled around with the lon values, but this is going to be a hassle
>> trying to eyeball on a plot what lat lon values to use to compensate for
>> what the script is doing
>>
>> Thank you again for all your help!
>>
>> Sho Kawazoe
>>
>> ;*****************************************************
>> ; Read original NetCDF file.
>> ;*****************************************************
>>
>> fin = addfile ("pr_MM5I_1979010103.nc","r")
>>
>> latS = 37
>> latN = 47
>> lonW = 269
>> lonE = 279
>>
>> lat = fin->lat ; read 2-D coordinates
>> lon = fin->lon
>> time = fin->time
>> Lambert_Conformal = fin->Lambert_Conformal
>>
>> printMinMax(lat, True)
>> printMinMax(lon, True)
>>
>> ji = region_ind (lat, lon, latS, latN, lonW, lonE)
>> j1 = ji(0) ; lat start ; get subset indices
>> j2 = ji(1) ; lat last
>> i1 = ji(2) ; lon start
>> i2 = ji(3) ; lon last
>>
>> prc = fin->pr(:,i1:i2,j1:j2) ; read 3-D subset
>>
>> lats_region = lat(i1:i2,j1:j2) ; make matching subset
>> lons_region = lon(i1:i2,j1:j2) ; of 2-D coordinates
>>
>> ; printVarSummary (prc) ; check results before proceeding
>> ; printMinMax(prc,True)
>> print ("")
>> print ("Lat should be 37-47")
>> print ("Lon should be 269-279")
>> printMinMax(lats_region,True)
>> printMinMax(lons_region,True)
>>
>>
>> ;******************************************************************
>> ; Create a new nc output file
>> ;******************************************************************
>>
>> system("/bin/rm -f region_output.nc") ; remove any pre-existing
>> file
>> fout = addfile("region_output.nc","c")
>> filedimdef(fout,"time",-1,True) ; Makes time an UNLIMITED
>> dimension (as suggested by NCL)
>>
>> fout->pr = prc
>> fout->lats_region = lats_region
>> fout->lons_region = lons_region
>>
>> fout->time = time
>> fout->Lambert_Conformal = Lambert_Conformal
>>
>> end
>>
>>
>> Output:
>> (0) latitude: min=20.7095 max=69.9062 (min max of input lat)
>> (0)
>> (0) longitude: min=203.718 max=322.936 (min max of input lon)
>> (0)
>> (0) Lat should be 37-47
>> (0) Lon should be 269-279
>> (0)
>> (0) latitude: min=54.8557 max=64.8567 (min max of output lat)
>> (0)
>> (0) longitude: min=229.353 max=253.274 (min max of output lon)
>>
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************
; -------------- MAIN -----------------
fName = "/Users/shea/Data/netCDF/pr_MM5I_1979010103.nc"
latS = 37
latN = 47
lonW = 269
lonE = 279
f = addfile (fName, "r")
PRC = f->pr ; read entire for demo
LAT = f->lat ; lat(yc, xc) ; latitude: min=20.7095 max=69.9062
LON = f->lon ; lon(yc, xc) ; longitude: min=203.718 max=322.936
printMinMax(LAT,True) ; entire grid
printMinMax(LON,True)
nlml = region_ind(LAT, LON, latS, latN, lonW, lonE)
nlStrt = nlml(0) ; lat start
nlLast = nlml(1) ; lat last
mlStrt = nlml(2) ; lon start
mlLast = nlml(3) ; lon last
prc = f->pr(:,nlStrt:nlLast,mlStrt:mlLast) ; read 3-D subset
latr = LAT(nlStrt:nlLast,mlStrt:mlLast) ; region subset
lonr = LON(nlStrt:nlLast,mlStrt:mlLast) ; of 2-D coordinates
print ("")
print ("Lat should be approximately 37-47")
print ("Lon should be approximately 269-279")
printMinMax(latr,True)
printMinMax(lonr,True)
;************************************************
; Sizes of arrays [grids]
;************************************************
dimPRC = dimsizes(PRC)
ntim = dimPRC(0)
NLAT = dimPRC(1)
MLON = dimPRC(2)
dimprc = dimsizes(prc)
nlat = dimprc(1)
mlon = dimprc(2)
;************************************************
; For convenience .... convert to kg/(m2-s) to mm/day
;************************************************
PRC = PRC*86400
prc = prc*86400
PRC@units = "mm/day"
prc@units = "mm/day"
;************************************************
; create plots
;************************************************
wks = gsn_open_wks("ps" ,"MM5") ; ps,pdf,x11,ncgm,eps
;gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
colors = (/"white","black","Snow" \ ; "WhiteSmoke" \
,"PaleTurquoise","PaleGreen","SeaGreen3" ,"Yellow" \
,"Orange","HotPink","Red","Violet", "Purple", "Brown"/)
gsn_define_colormap(wks, colors) ; generate new color map
res = True ; plot mods desired
res@gsnMaximize = True ; uncomment to maximize size
res@gsnSpreadColors = True ; use full range of colormap
res@cnFillOn = True ; color plot desired
res@cnFillMode = "CellFill" ; Cell Mode
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False ; turn off contour labels
;res@lbLabelAutoStride = True ; NCL should choose the spacing
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/0.1,1,2.5,5,10,15,20,25,50,75/) ; "mm/day"
;************************************************
; associate the 2-dimensional coordinates to the variable for plotting
;************************************************
PRC@lat2d = LAT ; direct assignment
PRC@lon2d = LON
prc@lat2d = latr ; direct assignment
prc@lon2d = lonr
;************************************************
; Turn on lat / lon labeling ; projection
; char Lambert_Conformal ;
; Lambert_Conformal:grid_mapping_name = "lambert_conformal_conic" ;
; Lambert_Conformal:standard_parallel = 30., 60. ;
; Lambert_Conformal:longitude_of_central_meridian = -97. ;
; Lambert_Conformal:latitude_of_projection_origin = 47.5 ;
;************************************************
res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks
lc = f->Lambert_Conformal
res@mpProjection = "LambertConformal"
res@mpLambertParallel1F = lc@standard_parallel(0)
res@mpLambertParallel2F = lc@standard_parallel(1)
res@mpLambertMeridianF = lc@longitude_of_central_meridian
;************************************************
; Demo: one arbitrarily chosen time
;************************************************
dimpr = dimsizes(prc) ; dimensions of prc
ntim = dimpr(0) ; number of time steps
nt = 0
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = LAT(0,0)
res@mpLeftCornerLonF = LON(0,0)
res@mpRightCornerLatF = LAT(NLAT-1,MLON-1)
res@mpRightCornerLonF = LON(NLAT-1,MLON-1)
plot = gsn_csm_contour_map(wks,PRC(nt,:,:),res)
res@mpLeftCornerLatF = latr(0,0)
res@mpLeftCornerLonF = lonr(0,0)
res@mpRightCornerLatF = latr(nlat-1,mlon-1)
res@mpRightCornerLonF = lonr(nlat-1,mlon-1)
plot = gsn_csm_contour_map(wks,prc(nt,:,:),res)
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 12 08:42:29 2010
This archive was generated by hypermail 2.1.8 : Fri Mar 12 2010 - 09:11:56 MST
