Re: how to link data with its geographic location information

Hi Ming,
It is a bit hard to diagnose the problems without access to your data. I can make a few suggestions though concerning your script.
Rather than use a loop to convert your data from 1D to 2D it would be much more efficient to use the NCL function onedtond (http://www.ncl.ucar.edu/Document/Functions/Built-in/onedtond.shtml)
Another thing is to make sure that the corner lat/lon values are correct. If the data is anything like the NSIDC data there may be some fill values around the edges in the lat/lon grids -- possibly the cause
of your trouble. It won't be exactly the same, but you might want to take a look at the example for snow and ice data at http://www.ncl.ucar.edu/Applications/nic_ims.shtml

If these suggestions don't help, if you can upload the data and coordinate files to the ftp site (see http://www.ncl.ucar.edu/report_bug.shtml for how to do this), we can take a look and probably help you.
 -dave

On Oct 30, 2013, at 11:04 AM, Ming Chen <chenming@ucar.edu> wrote:

> Greetings,
>
> I am trying to plot snow cover data generated by Rutgers University
> Climate Lab. There are two data sets in ascii format. One is snow cover
> data, which consists of 2D array "snow(i,j)", where i=1, nlon, j=1,nlat.
> The other data set gives the (lat,lon) information for each grid box
> corresponding to snow(i,j). The data is polar sterographic projection
> of northern hemisphere with central meridian -80.
>
> I read the two data sets, rearrange the data to be "snow(nlat,nlon)",
> and gives the corresponding geographic location "lat(nlat,nlon)" and
> "lon(nlat,nlon)". But I am not sure how to link "snow(nlat,nlon)" to
> "lat(nlat,nlon)" and "lon(nlat,nlon)".
>
> When I run the script, the error message is as follows:
>
> Warning:MapSetTrans: map limits invalid - using maximal area
>
> The figure generated is wrong.
>
> The script is attached. Any help and suggestions are highly appreciated.
>
> Ming
>
> ; Script to calculate domain-average vertical profile of QVAPOR,
> ; Then plot vertical- time series of these profiles
> ;
>
> 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/wrf/WRFUserARW.ncl"
>
> begin
> ;
> nlon=89
> nlat=89
> nmon=12
> nyear=1
> data=new((/nlat,nlon,1,12/),float)
> lat2d=new((/nlat,nlon/),float)
> lon2d=new((/nlat,nlon/),float)
>
> type = "x11"
> wks = gsn_open_wks(type,"polar")
> gsn_define_colormap(wks,"wh-bl-gr-ye-re")
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; read lat-lon info
> flnm = "rucl-gridboxes-center.txt"
> point = asciiread(flnm, -1, "string")
> delim = " "
> nfields = str_fields_count(point(0), delim) ; Count the fields
> separated
> ; by one or more spaces.
> print(nfields) ; nfields = 6
>
> field = 1
> subStrings = str_get_field(point, field, delim)
> ii = stringtointeger(subStrings)
> ;
> subStrings = str_get_field(point, field, delim)
> jj = stringtointeger(subStrings)
> ;
> field = 4
> subStrings = str_get_field(point, field, delim)
> lat= stringtofloat(subStrings)
> ;
> field = 5
> subStrings = str_get_field(point, field, delim)
> lon= stringtofloat(subStrings)
> ;
> field = 6
> subStrings = str_get_field(point, field, delim)
> ;lon = where(subStrings.eq."W", -lon, lon)
> lon = where(subStrings.eq."W", 360-lon, lon)
> ;
> ; rearrange the lat-lon to make it to 2D array
> ntot=dimsizes(lon)
> do i=0,ntot-1
> ix=mod(ii(i),90)-1
> do j=0,ntot-1
> jy=jj(j)-1
> lat2d(jy,ix)=lat(j)
> lon2d(jy,ix)=lon(j)
> end do
> end do
> print ( " REARRANGE END" )
>
> ; read snow data
> CASEDir = "/glade/scratch/chenming/WRFHELP/WQG/"
> FILES = systemfunc (" ls -1 " + CASEDir + "monmtx-1966-201231.cdr.snw")
> time = asciiread(FILES, -1, "string")
> nrow = dimsizes(time)
> print (nrow )
> ncount=0
> do irow=0,nrow-1
> if(mod(irow, 90).eq.0) then
> year=stringtointeger(str_get_cols(time(irow), 0,3))-1966
> mon=stringtointeger(str_get_cols(time(irow), 4,5))-1
> print (year )
> print (mon )
> ncount=ncount+1
> else
> newrow=irow-1-(ncount-1)*90
> print ( irow + " count= "+ ncount + " new row= " + newrow )
> do ilon=0,nlon-1
> test=stringtofloat(str_get_cols(time(irow), ilon*3,ilon*3+2))
> data(newrow, ilon,year,mon)=test ; data need to be
> array(nlat, nlon)
> end do
> end if
> end do
>
> ; now we need to define the lat & lon for data
> ; data!0 = "lat" ; This is to define the dimesnion names of the data
> ; data!1 = "lon"
> ; data&lat = lat2d ; This is to give the latitude of the data ( lat
> must be a 1-D array holding the latitude!)
> ; data&lon = lon2d
> ; lat@unit = " degrees_north"
> ; lon@unit = " degrees_east"
> ; print (data)
>
> ;start plotting
> res = True
> res@gsnAddCyclic = False ; regional data: not cyclic
> res@tfDoNDCOverlay = True ; set True for native projection
> res@gsnPolar = "NH" ; specify the hemisphere
> res@mpProjection = "Stereographic"
> res@mpLimitMode = "Corners"
> res@mpLeftCornerLatF = lat2d(0,0)
> res@mpLeftCornerLonF = lon2d(0,0)
> res@mpRightCornerLatF = lat2d(nlat-1,nlon-1)
> res@mpRightCornerLonF = lon2d(nlat-1,nlon-1)
> res@mpCenterLonF = -80
> res@mpCenterLatF = 90
>
> res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last
> res@mpFillDrawOrder = "PreDraw"
> res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
> res@mpFillOn = False ; turn off map fill
>
> res@gsnSpreadColors = True ; use full range of colormap
> res@cnFillOn = True ; color plot desired
> res@cnLinesOn = False ; turn off contour lines
> res@cnLineLabelsOn = False ; turn off contour labels
>
> ; plot = gsn_csm_contour_map_polar(wks,data(:,:,0,11),res) ;
> create the plot
> plot = gsn_csm_contour_map(wks,data(:,:,0,11),res) ; create the plot
>
> end
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Oct 30 12:57:18 2013