Re: reading ascii polar stereo and plotting?

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Wed, 25 Jun 2008 17:30:36 -0600

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