Hello
I am having trouble regridding MODIS data onto a 0.5 degree grid. I am using the MOD16 ET data, NCL version 6.1.2. I have generated a weighting file with ESMF_regrid_gen_weights, but when I attempt to interpolate, all of my data is missing.
I am using curvilinear_to_SCRIP for the src Grid. The MOD16 data is in a (1200x1200) array, from a 1 km resolution Swath. Then I use rectilinear_to_SCRIP to generate the destination grid, which I would like to be 0.5 degree over the region of the swath. Then I use ESMF_regrid_gen_weights and ESMF_regrid_with_weights to interpolate. I am not using ESMF_regrid right now because I wanted to try doing this step-by-step, to see where the problem occurs. I put some print statements within the procedure ESMF_regrid_with_weights, and my srcData is okay, but the dstData is all missing. This seems to occur during the function sparse_matrix_mult. One thing I've noticed is that some of my weights are negative - is this okay? Otherwise I guess I have done something wrong with treating missing data (a lot of the original data is missing - of the 1200x1200 array, there are only 66 non-missing values). Any guidance you can give is greatly appreciated!
I will include my code, and the print statements when I run it (in between the "--snip--" lines). I will also upload the following to the ftp (the source SCRIP file and weight file are very large though, so I'm having trouble uploading them):
MOD16A2.A2000M07.h00v09.105.2013119144046.hdf
MOD16A2.A2000M07.h00v09.105.2013119144046.he2
PET0.RegridWeightGen.Log
Rect.nc (the destination grid SCRIP file)
modis_regrid_esmf.ncl
Thanks
Anna
load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;shortened
begin
;SETTINGS
;---Output (and input) files
srcGridName = "MOD16_SCRIP.nc"
dstGridName = "Rect.nc"
wgtFile = "MOD_to_Rect.nc"
;*********************************
; Read in MODIS HDF files
;*********************************
hdf_file = addfile("MOD16A2.A2000M07.h00v09.105.2013119144046.hdf","r")
data_info=hdf_file->ET_1km
; Print information about the specific dataset to check scale and offset.
printVarSummary(data_info)
--snip--
Variable: data_info
Type: short
Total Size: 2880000 bytes
1440000 values
Number of Dimensions: 2
Dimensions and sizes: [YDim_MOD_Grid_MOD16A2 | 1200] x [XDim_MOD_Grid_MOD16A2 | 1200]
Coordinates:
Number Of Attributes: 10
scale_factor : 0.1
scale_factor_err : 0
add_offset : 0
add_offset_err : 0
calibrated_nt : 22
valid_range : ( -32767, 32700 )
_FillValue : 32767
long_name : MOD16A2 -- MODIS Gridded 1KM monthly Composite Evapotranspiration (ET)
units : kg/m^2/mon
hdf_name : ET_1km
--snip--
;READ THE ACTUAL DATA
eos_file = addfile("MOD16A2.A2000M07.h00v09.105.2013119144046.he2","r")
data_unscaled=eos_file->ET_1km_MOD_Grid_MOD16A2
data_unscaled@_FillValue = data_info@_FillValue
; Apply offset and scale using the formula in
; http://modis-atmos.gsfc.nasa.gov/MOD07_L2/format.html.
data_valid = where(data_unscaled.gt.data_info@valid_range(0) .and. \
data_unscaled.lt.data_info@valid_range(1), \
data_unscaled, data_info@_FillValue);
data = data_info@scale_factor * (data_valid - data_info@add_offset)
printMinMax(data,True)
printVarSummary(data)
--snip--
(0) min=75.8 max=143.3
Variable: data
Type: double
Total Size: 11520000 bytes
1440000 values
Number of Dimensions: 2
Dimensions and sizes: [1200] x [1200]
Coordinates:
Number Of Attributes: 1
_FillValue : 32767
--snip--
lat2d = eos_file->GridLat_MOD_Grid_MOD16A2
lon2d = eos_file->GridLon_MOD_Grid_MOD16A2
lon2d = where(lon2d.lt.0,360+lon2d,lon2d)
lat2d(:,:) = lat2d(::-1,::-1) ;reorder so lat is S->N
if (any(data.ne.data@_FillValue)) then
;----------------------------------------------------------------------
; Step 1, part 1
; Write MOD16 grid description to SCRIP file
;----------------------------------------------------------------------
Opt = True
Opt@ForceOverwrite = True
Opt@Title = "MOD16 grid"
;if the following is uncommented, I get no weights in the weight file
;from regrid_gen_weights
;Opt@Mask2D = where(.not.ismissing(data),1,0)
curvilinear_to_SCRIP(srcGridName,lat2d,lon2d,Opt)
;---Clean up
delete(Opt)
;----------------------------------------------------------------------
; Step 1, part 2
; Write rectilinear grid description to SCRIP file
;----------------------------------------------------------------------
;SET UP THE LAT/LON OF THE DESTINATION GRID
latS = doubletoint(min(lat2d))
latN = doubletoint(max(lat2d))
lonW = doubletoint(min(lon2d))
lonE = doubletoint(max(lon2d))
if(latS.gt.min(lat2d)) then
latS=latS-1
end if
if(latN.lt.max(lat2d)) then
latN=latN+1
end if
if(lonW.gt.min(lon2d)) then
lonW=lonW-1
end if
if(lonE.lt.max(lon2d)) then
lonE=lonE+1
end if
if(latN.lt.0 .and. latS.lt.0) then
nlat = (abs(latS)-abs(latN))*2 + 1
else if (latS.lt.0) then
nlat = (abs(latS)+abs(latN))*2 + 1
end if
end if
nlon = (lonE-lonW)*2 + 1
lat1d = fspan(latS, latN,nlat)
lon1d = fspan(lonW,lonE,nlon)
dims = dimsizes(data)
print("Original range of lat/lon")
printMinMax(lat2d,True)
printMinMax(lon2d,True)
print("New range: ")
print(nlat+" "+latS+" "+latN+" "+nlon+" "+lonW+" "+lonE)
--snip--
(0) Original range of lat/lon
(0)
(0) latitude: min=-9.99583 max=-0.00416667
(0)
(0) longitude: min=177.23 max=189.996
(0) New range:
(0) 21 -10 0 27 177 190
--snip--
Opt = True
Opt@ForceOverwrite = True
Opt@LLCorner = (/latS,lonW/)
Opt@URCorner = (/latN,lonE/)
;Opt@Debug = True
rectilinear_to_SCRIP(dstGridName,lat1d,lon1d,Opt)
;---Clean up
delete(Opt)
;----------------------------------------------------------------------
; Step 2
; Generate interpolation weights for WRF grid to Rectilinear grid
;----------------------------------------------------------------------
print("Generate interpolation weights for MOD grid to Rect grid")
Opt = True
Opt@InterpMethod = "Patch"
Opt@SrcRegional = True
Opt@DstRegional = True
Opt@PrintTimings = True
Opt@ForceOverwrite = True
Opt@Debug = True
ESMF_regrid_gen_weights(srcGridName, dstGridName, wgtFile, Opt)
--snip--
(0) Generate interpolation weights for MOD grid to Rect grid
(0) ESMF_regrid_gen_weights: number of processors used: 1
(0) --------------------------------------------------
(0) ESMF_regrid_gen_weights: the following command is about to be executed on the system:
(0) 'ESMF_RegridWeightGen --source MOD16_SCRIP.nc --destination Rect.nc --weight MOD_to_Rect.nc --method patch --src_regional --dst_regional -i --64bit_offset'
(0) --------------------------------------------------
(0) ESMF_regrid_gen_weights: output from 'ESMF_RegridWeightGen':
(0) Starting weight generation with these inputs:
(1) Source File: MOD16_SCRIP.nc
(2) Destination File: Rect.nc
(3) Weight File: MOD_to_Rect.nc
(4) Source File is in SCRIP format
(5) Source Grid is a regional grid
(6) Source Grid is a logically rectangular grid
(7) Destination File is in SCRIP format
(8) Destination Grid is a regional grid
(9) Destination Grid is a logically rectangular grid
(10) Regrid Method: patch
(11) Pole option: NONE
(12) Ignore unmapped destination points
(13) Output weight file in 64bit offset NetCDF file format
(14)
(15) Completed weight generation successfully.
(16)
(0) ESMF_regrid_gen_weights: 'ESMF_RegridWeightGen' was successful.
(0) --------------------------------------------------
--snip--
;---Clean up
delete(Opt)
;----------------------------------------------------------------------
; Step 3
; Interpolate data from WRF to regional rectilinear grid.
; Regrid two datasets: TMP and U_GRD.
;----------------------------------------------------------------------
;---Interpolate temperature and u
Opt = True
Opt@Debug = True
Opt@DstGridType = "rectilinear"
data_regrid = ESMF_regrid_with_weights(data,wgtFile,Opt)
--snip--
(0) ESMF_regrid_with_weights: regridding using interpolation weights ...
(0) ESMF_regrid_with_weights: warning: destination grid is not
(0) completely covered by the source grid. This is not an error.
(0) It just means your destination grid covers a larger area
(0) than your source grid.
(0) ESMF_regrid_with_weights: Source Grid:
(0) rank: 2
(0) dimensions: 1200 1200
(0) original source rank: 2
(0) latitude min/max: -9.99583/-0.00416667
(0) longitude min/max:177.23/189.996
(0) ESMF_regrid_with_weights: Destination Grid:
(0) dimensions: 21 27
(0) latitude min/max: -10/0
(0) longitude min/max:177/190
(0) ESMF_regrid_with_weights: retrieving interpolation weights ...
(0) ESMF_regrid_with_weights: calling sparse_matrix_mult to apply weights...
(0) ESMF_regrid_with_weights: dstData
(0) Dimensions: 21 27
(0) minSrcData: 75.8
(0) maxSrcData: 143.3
(0) minDstData: 32767
(0) maxDstData: 32767
--snip--
;---Reset 0 values to missing values.
data_regrid@_FillValue = data_info@_FillValue
data_regrid = where(data_regrid.eq.0.0,data_regrid@_FillValue,\
data_regrid)
printVarSummary(data_regrid)
print("Check min/max of regridded and original data, respectively:")
printMinMax(data_regrid,True)
printMinMax(data,True)
--snip--
Variable: data_regrid
Type: double
Total Size: 4536 bytes
567 values
Number of Dimensions: 2
Dimensions and sizes: [lat | 21] x [lon | 27]
Coordinates:
lat: [ -10.. 0]
lon: [ 177.. 190]
Number Of Attributes: 3
_FillValue : 32767
remap : remapped via ESMF_regrid_with_weights: Bilinear remapping
missing_value : 32767
(0) Check min/max of regridded and original data, respectively:
(0)
(0) min=32767 max=32767
(0)
(0) min=75.8 max=143.3
--snip--
end if
end
Anna Harper, PhD
College of Engineering, Mathematics, and Physical Sciences
University of Exeter
Harrison Building, North Park Road
Exeter, EX4 4QF
Phone: +44 1392 725910
http://emps.exeter.ac.uk/mathematics/staff/ah431
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jun 5 16:01:27 2013
This archive was generated by hypermail 2.1.8 : Tue Jun 11 2013 - 12:03:58 MDT