Anji
Attached is a script that will do the plotting. It uses explicit
contour levels. For precipitation fields it is best to use a non-evenly
spaced scale. You may wish to change. I used some from
an existing script I had.
The key is to get rid of the Nans.
In the attached NCL script, I just created a temporary file using the
unix "sed" tool.
This changed all the NaN to -999.
[snip]
; change all NaN to -999; create local work file "FOO"
system("sed s/NaN/-999/g "+diri+fili +">| FOO")
data = asciiread("FOO",(/nRow,nCol/),"float")
lat2d = onedtond(data(:,0), (/nlat,mlon/))
lon2d = onedtond(data(:,1), (/nlat,mlon/))
prate = onedtond(data(:,4), (/nlat,mlon/))
;err = onedtond(data(:,5), (/nlat,mlon/))
delete(data) ; no longer needed
system("/bin/rm -f FOO") ; delete local temp file
Good luck
D
Anji Seth wrote:
> Hello NCL people,
>
> I have an ascii dataset of estimated climatological snow accumulation
> for Antarctica
> in polar stereographic coordinates and I would like to read the file
> and plot it up in NCL.
>
> I have 2 questions:
> 1. The NCL example for reading gridded ascii data appears to be for a
> regular grid.
> Does anyone have an example of how to read polar stereo data?
>
> The data looks like this:
>
> (1) Latitude (degrees)
> (2) Longitude (degrees)
> (3) X(m) - Polar Stereographic (latitude of true scale 71S, ellipsoid
> WGS84)
> (4) Y(m) - Polar Stereographic (latitude of true scale 71S, ellipsoid
> WGS84)
> (5) Accumulation Rate (kg/m2/a)
> (6) Estimated r.m.s. error at 100 km scale, as a percentage of column
> 5, (%)
> -39.36 317.77 -3949233 4350420 NaN NaN
> -39.51 317.60 -3949231 4325353 NaN NaN
> -39.65 317.44 -3949238 4300275 NaN NaN
> .
> .
> .
>
> 2. Once I have read the data into NCL can I use
> gsn_csm_contour_map_polar ?
> The resolution of the data is effectively 100 km. Would I have to
> regrid it to T85 before
> I could plot or is there a way to plot on the native grid?
>
> Suggestions would be much appreciated!
> Thanks very much. -Anji
>
>
>
>
> Anji Seth, Geography, University of Connecticut
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-- ====================================================== Dennis J. Shea tel: 303-497-1361 | P.O. Box 3000 fax: 303-497-1333 | Climate Analysis Section | Climate & Global Dynamics Div. | National Center for Atmospheric Research | Boulder, CO 80307 | USA email: shea 'at' ucar.edu | ======================================================
; Anji Seth: 06/24/2008 09:52 AM
; (1) Latitude (degrees)
; (2) Longitude (degrees)
; (3) X(m) - Polar Stereographic (latitude of true scale 71S,
; ellipsoid WGS84)
; (4) Y(m) - Polar Stereographic (latitude of true scale 71S,
; ellipsoid WGS84)
; (5) Accumulation Rate (kg/m2/a)
; (6) Estimated r.m.s. error at 100 km scale,
; as a percentage of column 5, (%)
; Sample
; -39.36 317.77 -3949233 4350420 NaN NaN
; -39.51 317.60 -3949231 4325353 NaN NaN
; -39.65 317.44 -3949238 4300275 NaN NaN
;============================================================
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"
begin
diri = "/project/cas/shea/ANJI_SETH/"
fili = "amsr_accumulation_map.dat"
nCol = 6
nRow = numAsciiRow(diri+fili)
; change all NaN to -999; create local work file "FOO"
system("sed s/NaN/-999/g "+diri+fili +">| FOO")
data = asciiread("FOO",(/nRow,nCol/),"float")
nlat = 158
mlon = 664
lat2d = onedtond(data(:,0), (/nlat,mlon/))
lon2d = onedtond(data(:,1), (/nlat,mlon/))
prate = onedtond(data(:,4), (/nlat,mlon/))
;err = onedtond(data(:,5), (/nlat,mlon/))
delete(data) ; no longer needed
system("/bin/rm -f FOO")
lat2d_at_long_name = "latitude"
lat2d_at_units = "degrees_north"
lon2d_at_long_name = "longitude"
lon2d_at_units = "degrees_east"
prate_at_long_name = "Accumulation Rate"
prate_at_units = "kg/m2/a"
prate@_FillValue = -999.
; associate 2d coordinates with the variable
prate_at_lat2d = lat2d
prate_at_lon2d = lon2d
printVarSummary(prate)
printMinMax( prate, True)
printMinMax( lat2d, True)
printMinMax( lon2d, True)
;************************************************
; create plot
;************************************************
pltType = "x11" ; ps, x11, eps, pdf
pltTitle = "amsr_prate"
wks = gsn_open_wks(pltType , pltTitle)
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
res = True ; plot mods desired
res_at_gsnPolar = "SH" ; specify the hemisphere
res_at_gsnMaximize = True ; make large id "ps" or "eps"
res_at_gsnPolarLabelSpacing = 90 ; how frequently to label
res_at_gsnSpreadColors = True
res_at_mpMaxLatF = -60.
res_at_cnFillOn = True ; color
res_at_cnLinesOn = False ; no contour lines
res_at_cnLineLabelsOn = False
res_at_mpFillOn = False
res_at_cnLevelSelectionMode = "ExplicitLevels"
res_at_cnLevels = (/10,25,50,75,100,150,200,300,500,1000/) ; mm
plot = gsn_csm_contour_map_polar(wks,prate,res) ; create the plot
; Trick: Don't plot anything less than 1
; not needed here
; prate = where(prate.le.1, prate@_FillValue, prate)
; plot = gsn_csm_contour_map_polar(wks,prate,res) ; create the plot
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jun 25 2008 - 17:30:36 MDT
This archive was generated by hypermail 2.2.0 : Thu Jun 26 2008 - 16:29:44 MDT