regridding mod10 data

Hi NCL-users,

I am having problems regridding MODIS swath data products to a regular lat/lon grid using rcm2rgrid. I use this function all the time with WRF data but cannot get it going for MODIS swath .hdf data. NCL does not crash when trying to evaluate rcm2rgrid command, it just hangs.

I have reduced the regrid area to a small region to minimise computational efforts and it still hangs. I have had problems with getting the geo lat/lon array in the correct order for rcm2rgrid previously when working in the southern hemisphere, so I suspect this is the source of the problem, however I believe I have conformed my data to the function requirements. A clean version of the script is below, any help would be well appreciated.

Kat

;name:
; kbormann_mod10.ncl
;
;purpose:
; this script reads mod10 fractional snow cover data from hdf files directly
; regrids to regular lat/lon grid and plots
;
;pre processing:
; need to get 5km swath lat/lons at 500m (hdf file only has 5km)
; you can use eos2dump program (hdfeos.org) to do a filedump of these values
; $cd /c/z3273429/bin/eos2dump/
; $./eos2dump_V1.0 -a2 /**pathtofile*/MOD10_L2.A2000217.0025.005.2007189173800.hdf > lon_MOD10_L2.A2000217.0025.005.2007189173800.hdf_output.txt
; $./eos2dump_V1.0 -a1 /**pathtofile*/MOD10_L2.A2000217.0025.005.2007189173800.hdf > lat_MOD10_L2.A2000217.0025.005.2007189173800.hdf_output.txt
;
;author:
; KBormann @ UNSW
;
;last modified:
; 20/11/2011
; 28/11/2011 - still debugging
;-

begin
print("____________________________________")

;---libraries
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

;---mod10 fractional snow cover product
  ;read data
  f = addfile("MOD10_L2.A2000217.0025.005.2007189173800.hdf","r")
  fsc_byt = f->Fractional_Snow_Cover
  lat_5km = f->Latitude
  lon_5km = f->Longitude
  print(" mod10 data read")

  ;process data format
  ;delete(fsc_byt@_FillValue)
  fsc = byte2flt(fsc_byt)
  fsc@_FillValue = 255
  printVarSummary(fsc)

;---mod10 geo coordinates
  ;currently lat/lon data at 5km resolution and fsc data at 500m resolution
  ;for information on how to obtain the lat/lon data, http://hdfeos.org/zoo/note_non_geographic.php
  ;follow example: http://www.hdfeos.org/zoo/NSIDC/MOD10_L2.A2000065.0040.005.2008235221207_snow_cover_P.ncl
  ;read output from eos2dump

  ;data sizes
  ds=dimsizes(fsc)
  nlon=ds(0) ;4060 (from datafield in hdf file)
  nlat=ds(1) ;2708 (from datafield in hdf file)

  ;read data
  flat = "lat_MOD10_L2.A2000217.0025.005.2007189173800.hdf_output.txt"
  flon = "lon_MOD10_L2.A2000217.0025.005.2007189173800.hdf_output.txt"
  lat2d = asciiread(flat,(/nlon,nlat/),"float")
  lon2d = asciiread(flon,(/nlon,nlat/),"float")
  print(" mod10 geo coordinates read")

  ;assign correct geo data to snow data
  fsc@lat2d=lat2d
  fsc@lon2d=lon2d

;---assign zoom region, don't need entire swath
   minlon = 148.2
   maxlon = 148.6
   minlat = -36.4
   maxlat = -36.0

;---plot swath data
  ;make maximum filesize larger
  setvalues NhlGetWorkspaceObjectId()
  "wsMaximumSize" : 200000000
  end setvalues

   wks = gsn_open_wks("x11", "swathplot")
   gsn_define_colormap(wks,"rainbow+white+gray") ; choose colormap

   res = True
   res@cnFillOn = True ; color fill
   res@cnFillMode = "CellFill" ; raster mode
   res@cnLinesOn = False ; turn off contour lines .
   res@cnLineLabelsOn = False ; turn off contour labels.
   res@gsnMaximize = True ; make plot large

   res@cnMissingValFillPattern= 2 ; missing value pattern is set to "SolidFill"
   res@cnMissingValFillColor = 239 ; white color for missing values
   res@cnFillColors = (/238,238,111,104,96,88,80,72,64,56,48,42,239,239,239,144,25,26,190,239/)
   res@tiMainString = "mod10 fsc swathplot"
   res@tiMainFontHeightF = 0.018
   res@tiXAxisString = "lon"
   res@tiYAxisString = "lat"
   res@tiXAxisFontHeightF = 0.015
   res@tiYAxisFontHeightF = 0.015
   res@gsnAddCyclic = False
   res@mpLimitMode = "LatLon"
   res@pmTickMarkDisplayMode = "Always"
   res@cnRasterSmoothingOn = False
   res@cnLevelSelectionMode = "ExplicitLevels" ; set explict contour levels
   res@cnLevels = (/0,10,20,30,40,50,60,70,80,90,100,200,201,211,225,237,239,250/)

   res@lbLabelPosition = "Center" ; label position
   res@lbLabelAlignment = "BoxCenters" ; label orientation
   res@lbLabelStrings = (/"-","0","10","20","30","40","50","60","70","80","90","100","200","201","211","225","237","239","250","+"/)
   res@lbLabelFontHeightF = 0.010
   res@lbOrientation = "Horizontal"
   res@gsnPaperOrientation = "Portrait"

   res@lbTitleString = (/"KEY: 0-100=snow fsc %, 200=missing data, 201=no decision,~C~211=night, 225=land, 237=inland water, 239=ocean, 250=cloud"/)
   res@lbTitlePosition = "Bottom"
   res@lbTitleFontHeightF = 0.011

   res@mpMinLonF = minlon
   res@mpMaxLonF = maxlon
   res@mpMinLatF = minlat
   res@mpMaxLatF = maxlat

   plot = gsn_csm_contour_map_ce(wks,fsc, res) ; create plot
   print(" swath data plotted")

;---regrid swath data to regular grid
   ;make regular grid coordinates roughly around zoom region
   ;TL = topleft, BR = bottomright
   TLlat = maxlat+0.5
   TLlon = minlon-0.5
   BRlat = minlat-0.5
   BRlon = maxlon+0.5
   print(" desired regrid bounds: "+TLlat+","+TLlon+","+BRlat+","+BRlon)

   ;generate rectilinear grid for subset region 500m res.
   deg500m = 0.004166667 ;500m expressed in decimal degrees
   nptsY500 = round(((TLlat - BRlat)/deg500m),3)
   nptsX500 = round(((BRlon - TLlon)/deg500m),3)
   rlat = fspan(BRlat, TLlat, nptsY500)
   rlon = fspan(TLlon, BRlon, nptsX500)
   print(" check rect grid coords are monotonically increasing...")
   print(" rlat(0)="+rlat(0)+", rlat(10)="+rlat(10))
   print(" rlon(0)="+rlon(0)+", rlon(10)="+rlon(10))

   ;assign attributes to new grid
   rlat@res = "500m intervals within subset"
   rlat@units = "degrees north"
   rlon@res = rlat@res
   rlon@units = "degrees east"

   ;check 2d coords are correctly defined
   lat2d=lat2d(::-1,::-1)
   lon2d=lon2d(:,::-1)

   ;need to adjust fsc to reflect these changes (not yet as regridding not working)

   ;rcm2rgrid help says south-to-north, require [ lat2d(j+1,0) > lat2d(j,0) ]
   print(" check 2d grid coords are south-to-north")
   print(" lat2d(0,0)="+lat2d(0,0)+" to lat2d(0,10)="+lat2d(0,10))
   print(" lat2d(0,0)="+lat2d(0,0)+" to lat2d(10,0)="+lat2d(10,0))
   ;east-to-west
   print(" check 2d grid coords are east-to-west")
   print(" lon2d(0,0)="+lon2d(0,0)+", lon2d(0,10)="+lon2d(0,10))
   print(" lon2d(0,0)="+lon2d(0,0)+", lon2d(10,0)="+lon2d(10,0))

   ;check desired zoom coordinates are within swath geolocation
   print(" check desired grid lats are actually in swath")
   if ( (min(lat2d).le.min(rlat)).and.(max(lat2d).ge.max(rlat)) ) then
      print(" Y")
   else
      print(" N")
     exit
   end if
   print(" check desired grid lons are actually in swath")
   if ( (min(lon2d).le.min(rlon)).and.(max(lon2d).ge.max(rlon)) ) then
      print(" Y")
   else
      print(" N")
     exit
   end if

   ;assign dimension names in case this is important to rcm2rgrid
   printVarSummary(lon2d)
   printVarSummary(lat2d)
   fsc!0 = "latitude"
   fsc!1 = "longitude"
   printVarSummary(fsc)

   ;regrid subset data to 500m rectilinear grid
   print(" regridding mod10 fsc data to 500m rectilinear grid at zoom...")
   rgridded = rcm2rgrid_Wrap(lat2d, lon2d, fsc, rlat, rlon,1)
   printVarSummary(rgridded)

;the whole point of this script is to get swath mod10 data to match landsat data 2d grid
;can't go directly from a 2d grid to another 2d grid in ncl. need to get to regular grid first.

;in the future regrid to exact landsat coordinates
      ;regrid subset data to landsat fsc coordinates
       ;print(" regridding mod10 fsc data to landsat fsc 2d grid...")
       ;fsc_lan_grid = rgrid2rcm_Wrap(lat2d, lon2d, rgridded, lat2d_lan, lon2d_lan,1)
      ;printVarSummary(fsc_lan_grid)

print("____________________________________")
end

Kat Bormann
PhD Research Candidate
Climate Change Research Centre
University of New South Wales
Kensington Campus, Sydney, NSW
Email: k.bormann@student.unsw.edu.au
http://web.science.unsw.edu.au/~katbormann/
http://hydrology.unsw.edu.au/student/bormann/
Mob: 0402 573 380

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 28 19:20:50 2011