[SNIP]
>
> I'm not familiar with if the valid_range attributes are actually supposed to be used to scale the data. I think they might just be informational. I've CC'ed Dennis as he may know more about this.
>
These are "informational". The should not be used to scale the data.
> Dennis, this is how the data is being scaled in the script. Note the first read is to get the proper _FillValue (I think). [I'm off to a meeting at 10, then in my office around 11:30 or so).
>
In my experience, there are two methods used to descale values.
If xs is the short number
(1) x = xs*scale_factor + add_offset
(2) x = (xs - add_offset)*scale_factor
Method 1 is required by the netCDF COARDS and CF Conventions.
NASA labs use 1 or 2. In most cases, there is no documentation
to users which to uses
f = addfile("MOD16A3.A2010365.h19v06.105.2011050092657.hdf","r")
data_info= f ->ET_1km
data = data_info@scale_factor * (data_info - data_info@add_offset)
data@long_name = data_info@long_name
data@units = data_info@units
data@_FillValue = 1e20 ; cange _FillValue
copy_VarCoords(data_info, data)
printVarSummary(data)
print("data: min="+min(data)+" max="+max(data))
>
> ;*********************************
> ; 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)
>
> ;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)
>
>
>
> On Jun 5, 2013, at 4:01 PM, "Harper, Anna" <A.Harper@exeter.ac.uk> wrote:
>
>> 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
>
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jun 6 15:13:03 2013
This archive was generated by hypermail 2.1.8 : Tue Jun 11 2013 - 12:03:58 MDT