Re: plotting CALIPSO data

Hi Prabhakar,

What about simply using "ind" to cull out the indices that contain the area of interest?

  minlat = 20 ; or some val
  maxlat = 50 ; or some val
  minhgt = 0
  maxhgt = 100 ;

  ii = ind(lat2d.ge.minlat.and.lat2d.le.maxlat.and.height.ge.minhgt.and.height.le.maxhgt)

and then subscript "lat2d", "height" and "var" with ii:

  res2@sfXArray = lat2d(ii)
  res2@sfYArray = height(ii)
  plot(1) = gsn_csm_contour(wks,var(ii), res2)

In order for this all to work, var, height, and lat2d must all be 1D arrays of the same length.
You might want to rename "lat2d" to "lat1d" so that you don't mistakenly treat it as a 2D array.

You mentioned using lat/lon to cull out the area of interest, but you are plotting lat/height. You can
replace the above "height" with your longitude array if needed.

--Mary

On May 2, 2010, at 1:00 PM, p s wrote:

> Hi
> I am getting into all sorts of trouble plotting the CALIPSO VFM data. The arrays are very large 4400*5000, and I get error messages:
> (0) 0 CAL_LID_L2_VFM-Prov-V2-01.2006-06-13T07-20-12ZD.hdf
> warning:ScalarFieldSetValues: coordinate array sfXArray requires 24964297 elements: defaulting
> warning:ScalarFieldSetValues: coordinate array sfYArray requires 24964218 elements: defaulting
>
> How do we just subset the plot into small region based on lat and long, so that I can plot it easily.
> The lat , long are in separate arrays and are 1-D, so I could not use closest_val to get the location of indices from where to subset as well.
>
> Looking forward for you kind suggestion.
> Regards,
> Prabhakar
>
> PS: I have also attached the code I am using.
>
> ;hdfreadCALIOP.ncl
> ;HERE WE READ AND PLOT THE VFM V2.0 products
>
> 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 = "/Data/CALIPSO/"
> fili = systemfunc("cd "+diri+" ;ls CAL* " )
> nfil = dimsizes( fili )
>
> do nf = 0,nfil-1 ;BIG LOOP
> print(nf+" "+fili(nf))
> f = addfile(diri+fili(nf), "r")
>
> var = f->Feature_Classification_Flags
>
> N = dimsizes(var)
> height = ispan(0,1000,1)
>
> ; GEOLOCATION ARRAYS
>
> lat2d = f->Latitude(fakeDim0|:,fakeDim1|0)
> lon2d = f->Longitude(fakeDim2|:,fakeDim3|0)
>
> ; PLOTTING
>
> if (nf.lt.10) then
> wks = gsn_open_wks ("eps", "vfm_00"+sprintf("%1.0f",nf))
> else
> if (nf.lt.100) then
> wks = gsn_open_wks ("eps", "vfm_0"+sprintf("%2.0f",nf))
> else
> wks = gsn_open_wks ("eps", "vfm_"+sprintf("%3.0f",nf))
> end if
> end if
>
> ; colors = (/"red","RoyalBlue","LightSkyBlue","darkorange","yellow","yellowgreen","thistle","black"/)
> ; gsn_define_colormap(wks,colors)
>
> ; START PANEL
>
> plot = new(2,graphic)
>
> res = True
>
> res@gsnDraw = False ; don't draw
> res@gsnFrame = False ; don't advance frame
> res@gsnMaximize = True ; maximize pot in frame
>
> ; SPECIFY SUBPLOT SPECIFIC GEOLOCATION AND TITLE
> res1 = res
> res1@tiMainString = "VFM"
>
> res1@mpOutlineBoundarySets = "National"
> res1@pmTickMarkDisplayMode = "Always"
> res1@mpDataBaseVersion = "MediumRes" ; Higher res coastline
> res1@mpLimitMode = "LatLon"
> res1@mpMinLatF = 20.
> res1@mpMaxLatF = 40.
> res1@mpMinLonF = 60.
> res1@mpMaxLonF = 90.
>
> ;
> res2 = res
> res2@tiMainString = "VFM"
> res2@cnFillOn = True ; color Fill
> res2@cnFillMode = "RasterFill" ; Raster Mode
> res2@cnRasterSmoothingOn = True
> ; res2@cnRasterMinCellSizeF = 0.0005
> res2@cnLinesOn = False ; Turn off contour lines
> res2@cnLineLabelsOn = False ; Turn off contour lines
>
> res2@sfXArray = lat2d
> res2@sfYArray = height
> res2@trGridType = "TriangularMesh"
> res2@gsnAddCyclic = False ; Data is not cyclic
> ;
> res2@lbLabelStrings = ispan(0,7,1) ; 0, 1, 2, 3, 4
> res2@lbLabelPosition = "Center" ; label position
> res2@lbLabelAlignment = "BoxCenters"
> res2@lbTitleString =" "; "0=low/no confidence, 1=clear air, 2=cloud , 3=aerosol, 4=stratospheric features,/
> ; 5=surface, 6=subsurface, 7=no signal"
> res2@lbTitlePosition = "Bottom"
> res2@lbTitleFontHeightF = 0.0125 ; default 0.025
>
> ;
> plot(0) = gsn_csm_map_ce(wks,res1) ; create map
> pres = True
> pres@gsLineThicknessF = 3.0
> pres@gsLineColor ="red"
> line = gsn_add_polyline(wks,plot(0),lon2d,lat2d,pres)
>
> plot(1) = gsn_csm_contour(wks,var, res2)
>
> resP = True
> resP@txString = date_str
> gsn_panel(wks,plot,(/2,1/),resP)
> draw(plot)
> frame(wks)
> end do
> 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 Mon May 3 10:52:06 2010