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
This archive was generated by hypermail 2.1.8 : Wed Nov 30 2011 - 19:52:47 MST