Problems of making plots with station data

> Hi,NCL users,
> I'd like to make plots with station data (ASII format) compare the
> output from RegCM3,the data is like:
> ID1 LAT1 LON1 PRE1
> ID2 LAT2 LON2 PRE2
> ...
> I have written the ncl script like the fllowing:
> You can find from the script that I used the model grid information to
> replace station information in stn_latlon.dat and for lambert
> projection since I haven't the stn_latlon.dat and don't know how to
> make it.In my option, that's the same thing.But oberviously, I can not
> make that in the correct way, the station mark can't correctly marked with
> the map,it seems like the rain and the station point has the different
> map or projection.One more thing, I want to mask the rain with the
> landuse value, but with the station data,it seems I can't do that.Could
> you give me some suggestions?
> Thanks a lot.
>
> Yan
>
> ---------------------------------------------------------------------------
> ; pathnames must be changed based on local installation of NCL
> load "/home/RAID2-D1/NCL/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "/home/RAID2-D1/NCL/lib/ncarg/nclscripts/csm/shea_util.ncl"
> load "/home/RAID2-D1/NCL/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "/home/RAID2-D1/NCL/lib/ncarg/nclscripts/csm/contributed.ncl"
>
> begin
> ; Set path to HEAD_OUT.nc
> diri1 = "/DATA/NCL/"
> fili1 = "out_head.nc"
> ; Pointers to input files
> f1 = addfile(diri1+fili1,"r")
> f2 = "rain.dat"
> ; Read in each row as a single string, and then convert to a character
> ; array so we can parse it.
> ;
> data = asciiread("rain.dat", -1, "string") ; -1 means read all rows.
> cdata = stringtochar(data)
> ;
> ; The first row is just header information, so we can discard this.
> ; That is, we can start with the second row, which is represented
> ; by index '1'.
> ;
> ; The latitude values fall in columns 6-12 (indices 7:13)
> ; The longitude values fall in columns 13-21 (indices 14:22)
> ; The pwv data values fall in columns 22-31 (indices 23:)
> ;
> lat = stringtofloat(charactertostring(cdata(1:,7:13)))
> lon = stringtofloat(charactertostring(cdata(1:,14:22)))
> rain = stringtofloat(charactertostring(cdata(1:,23:)))
> ;
>
> delete(data) ; Remove arrays we don't need anymore.
> delete(cdata)
> lat2d = f1->xlat(0,:,:) ; (time,lat,lon)
> lon2d = f1->xlong(0,:,:) ; (time,lat,lon)
> lu=f1->landuse(0,:,:);(time,lat,lon)
> ; Get dimensions from variable
> dimvar = dimsizes(var)
> jlat = dimvar(0)
> ilon = dimvar(1)
>
> wks = gsn_open_wks ("ps", "RegCM_rain") ; open workstation
> gsn_define_colormap(wks,"gui_default")
>
> res = True ; plot mods desired
> ; !!!!! any plot of data that is on a native grid, must use the "corners"
> ; method of zooming in on map.
>
> res_at_mpLimitMode = "Corners" ; choose range of map
> res_at_mpLeftCornerLatF = lat2d(0,0)
> res_at_mpLeftCornerLonF = lon2d(0,0)
> res_at_mpRightCornerLatF = lat2d(jlat-1,ilon-1)
> res_at_mpRightCornerLonF = lon2d(jlat-1,ilon-1)
>
> ; res_at_mpProjection = "LambertConformal"
> ; res_at_mpLambertParallel1F = tlat1
> ; res_at_mpLambertParallel2F = tlat2
> ; res_at_mpLambertMeridianF = clon
>
> ; res_at_mpLimitMode = "Corners"
> ; res_at_mpLeftCornerLatF = lat2d(20,18)
> ; res_at_mpLeftCornerLonF = lon2d(20,18)
> ; res_at_mpRightCornerLatF = lat2d(58,52)
> ; res_at_mpRightCornerLonF = lon2d(58,52)
> ; The following 4 pieces of information are REQUIRED to properly display
> ; data on a native lambert conformal grid. This data should be specified
> ; somewhere in the model itself.
>
> res_at_mpProjection = "LambertConformal"
> res_at_mpLambertParallel1F = 30.
> res_at_mpLambertParallel2F = 60.
> res_at_mpLambertMeridianF = 105
>
> ; usually, when data is placed onto a map, it is TRANSFORMED to the
> specified
> ; projection. Since this model is already on a native lambert conformal
> grid,
> ; we want to turn OFF the tranformation.
>
> res_at_tfDoNDCOverlay = True
> res_at_gsnFrame = False ; So we can draw markers
> res_at_gsnMaximize = True
>
> res_at_gsnSpreadColors = True
> res_at_gsnSpreadColorEnd = 90
> res_at_pmTickMarkDisplayMode = "Always"
>
> res_at_trGridType = "TriangularMesh" ; The default if
> you
> ; have 1D data
>
> res_at_cnLineLabelPlacementMode = "Constant"
> res_at_cnLineDashSegLenF = 0.3
>
> res_at_cnLevelSelectionMode = "ManualLevels"
> res_at_cnMinLevelValF = 15 ; 15.25
> res_at_cnMaxLevelValF = 50 ; 49.75
> res_at_cnLevelSpacingF = 0.25
>
> res_at_cnFillOn = True
> res_at_cnLinesOn = True
> res_at_cnLineLabelsOn = True
> res_at_cnLevelFlags = new(139,"string")
> res_at_cnLevelFlags(:) = "NoLine"
> res_at_cnLevelFlags(0::20) = "LineAndLabel"
>
> res_at_lbLabelStride = 20
> res_at_lbBoxLinesOn = False
> res_at_tiMainString = "GPS PWV (18Z)"
>
> res_at_sfXArray = lon
> res_at_sfYArray = lat
>
> res_at_mpFillOn = False
> res_at_mpOutlineDrawOrder = "PostDraw"
> res_at_mpFillDrawOrder = "PreDraw"
> res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"
> res_at_mpUSStateLineColor = "Gray10"
> res_at_mpUSStateLineDashPattern = 2
>
> ;
> ; Draw the map (frame won't get advanced because gsnFrame was set to
> False).
> ;
> map = gsn_csm_contour_map(wks,rain,res)
>
> ;
> ; Draw markers on the plot in the lat/lon locations.
> ;
> mkres = True
> mkres_at_gsMarkerIndex = 17 ; Filled circle
> mkres_at_gsMarkerSizeF = 0.03
>
> gsn_polymarker(wks,map,lon,lat,mkres)
>
> frame(wks) ; Now advance the frame.
>
>
> end
>
>
>
>
>
>

-- 
Yan Bao
Earth System Physics Section
The Abdus Salam International Centre for Theoretical Physics
P.O. BOX 586, (Strada Costiera 11 for courier mail)
34100 Trieste, Italy
Phone: + 39 040 2240246
Fax: + 39 040 2240449
email: ybao_at_ictp.it
_______________________________________________
ncarg-talk mailing list
ncarg-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncarg-talk