Re: regridding mod10 data

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Wed Nov 30 2011 - 08:32:10 MST

Kat and others,

We are hoping to add some new regridding software to NCL in the next release, using the ESMF regridding application (http://www.earthsystemmodeling.org/).

This work is the result of a successful student project by Mohammad Abouali, which we are generalizing for inclusion in NCL. We are also working on a suite of examples.

I have Kat's files, and will look into using this software to regrid the data.

Meanwhile, if anybody else has some data they'd like to regrid, please send me an email. It would help if you can provide a sample dataset, and a clear description (or a dataset) of the type of grid you want to interpolate to.

--Mary

On Nov 28, 2011, at 7:20 PM, Kat Bormann wrote:

> 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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Nov 30 08:32:20 2011

This archive was generated by hypermail 2.1.8 : Wed Nov 30 2011 - 19:52:47 MST