how to link data with its geographic location information

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
Received on Wed Oct 30 11:04:53 2013