Re: MODIS HDF

Thanks to Mary, Dennis, and David.

Mary, it seems now I don't get an error but
when I either use ncview on test.nc or try to plot
the data, it looks like garbage. It is not plotting correctly.
What am I doing wrong?

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/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/csm/shea_util.ncl"
begin

; 1km annual npp modis on 1200x1200 tiles

a=addfile("Y2000/MOD17A3.A2000366.h12v04.305.2009342145314.hdf.hdfeos","r")

lat2d=a->GridLat_MOD_Grid_MOD17A3
lon2d=a->GridLon_MOD_Grid_MOD17A3

npp=short2flt(a->Npp_1km_MOD_Grid_MOD17A3)

npp!0="nlat"
npp!1="nlon"
npp@coordinates="lat2d lon2d"

lat2d!0="nlat"
lat2d!1="nlon"
lon2d!0="nlat"
lon2d!1="nlon"

delete(npp@projection)
delete(lon2d@projection)
delete(lat2d@projection)

printVarSummary(npp)

system("rm test.nc")
out=addfile("test.nc","c")
out->npp=npp
out->lat2d=lat2d
out->lon2d=lon2d

npp=mask(npp,npp.lt.0.,False)

wks = gsn_open_wks("ps","read_npp")
setvalues wks
  "wkColorMap" : "rainbow+white"
end setvalues
plot=new(1,graphic)

res = True
res@cnLinesOn = False
res@cnFillOn = True
res@cnLineLabelsOn = False

res@cnRasterModeOn=True

res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat2d(0,0)
res@mpLeftCornerLonF = lon2d(0,0)
res@mpRightCornerLatF = lat2d(1199,1199)
res@mpRightCornerLonF = lon2d(1199,1199)

res@mpProjection = "LambertConformal"
res@mpLambertParallel1F = 35.
res@mpLambertParallel2F = 45.
res@mpLambertMeridianF = -85.

res@tfDoNDCOverlay = True

res@mpGeophysicalLineColor = 2
res@mpPerimOn = True
res@mpGridLineDashPattern = 2
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpUSStateLineColor = "black"
res@pmTickMarkDisplayMode = "Always"

res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/1000.,2000.,3000.,4000.,5000.,6000.,7000.,8000./)
res@cnFillColors=(/237,200,187,177,238,97,80,47,20/)
res@tiMainString="NPP"
plot(0) = gsn_csm_contour_map(wks,npp,res)

end

On Aug 30, 2011, at 2:58 PM, Dennis Shea wrote:

> OK ...
>
> Also, Mike ...
>
> You do not need to include projection. Once you have the values
> and the lat/lon locations, you can use whatever
> projection you want.
>
> D
>
> On 8/30/11 1:55 PM, Mary Haley wrote:
>> The problem seems to be related to the fact that the attribute "projection" is (null). If I delete the projection attribute or set it to some string before I write it to the file, then I don't get a seg fault.
>>
>> delete(npp@projection)
>> delete(lon2d@projection)
>> delete(lat2d@projection)
>>
>> I've submitted a JIRA ticket, NCL-1190.
>>
>> --Mary
>>
>> On Aug 30, 2011, at 1:47 PM, Michael Notaro wrote:
>>
>>> Thank you.
>>>
>>>
>>> On Aug 30, 2011, at 2:40 PM, Mary Haley wrote:
>>>
>>>> I'll take a look at this and try to narrow it down for Dave.
>>>>
>>>> --Mary
>>>>
>>>> On Aug 30, 2011, at 1:27 PM, Dennis Shea wrote:
>>>>
>>>>> DB .. I put the script and file at:
>>>>>
>>>>> ftp ftp.cgd.ucar.edu
>>>>>
>>>>> cd pub/shea/NOTARO
>>>>> prompt
>>>>> mget *
>>>>> quit
>>>>>
>>>>>
>>>>> this will get a seg fault
>>>>>
>>>>>
>>>>> -------- Original Message --------
>>>>> Subject: Re: MODIS HDF
>>>>> Date: Tue, 30 Aug 2011 13:25:15 -0500
>>>>> From: Michael Notaro<mnotaro@wisc.edu>
>>>>> To: Dennis Shea<shea@ucar.edu>
>>>>> CC: ncl-talk NCL<ncl-talk@ucar.edu>
>>>>>
>>>>> Thanks Dennis and David.
>>>>>
>>>>> I am still unfortunately getting a segmentation fault with
>>>>> the edited code below.
>>>>> (I believe I correctly implemented your suggestions, but
>>>>> let me know if I messed up anywhere)
>>>>>
>>>>> By the way, how do I verify the projection type? I suspect
>>>>> it is sinusoidal but am not sure.
>>>>>
>>>>>
>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/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/csm/shea_util.ncl"
>>>>> begin
>>>>>
>>>>> ; 1km annual npp modis on 1200x1200 tiles
>>>>>
>>>>> a=addfile("/Users/notaro/downloads/Y2000/MOD17A3.A2000366.h11v11.305.2009342143541.hdf.hdfeos","r")
>>>>>
>>>>> lat2d=a->GridLat_MOD_Grid_MOD17A3
>>>>> lon2d=a->GridLon_MOD_Grid_MOD17A3
>>>>>
>>>>> npp=short2flt(a->Npp_1km_MOD_Grid_MOD17A3)
>>>>>
>>>>> npp!0="nlat"
>>>>> npp!1="nlon"
>>>>> npp@coordinates="lat2d lon2d"
>>>>> npp@projection="sinusoidal"
>>>>>
>>>>> lat2d!0="nlat"
>>>>> lat2d!1="nlon"
>>>>> lon2d!0="nlat"
>>>>> lon2d!1="nlon"
>>>>>
>>>>> printVarSummary(npp)
>>>>>
>>>>> system("rm test.nc")
>>>>> out=addfile("test.nc","c")
>>>>> out->npp=npp
>>>>> out->lat2d=lat2d
>>>>> out->lon2d=lon2d
>>>>>
>>>>> end
>>>>>
>>>>>
>>>>>
>>>>> Variable: npp
>>>>> Type: float
>>>>> Total Size: 5760000 bytes
>>>>> 1440000 values
>>>>> Number of Dimensions: 2
>>>>> Dimensions and sizes: [nlat | 1200] x [nlon | 1200]
>>>>> Coordinates:
>>>>> Number Of Attributes: 5
>>>>> _FillValue : 1e+20
>>>>> coordinates : lat2d lon2d
>>>>> hdfeos_name : Npp_1km
>>>>> projection : sinusoidal
>>>>> unsigned : True
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Aug 30, 2011, at 1:16 PM, Dennis Shea wrote:
>>>>>
>>>>>>
>>>>>>
>>>>>> On 8/30/11 12:14 PM, Dennis Shea wrote:
>>>>>>>
>>>>>>>
>>>>>>> npp!0="nlat"
>>>>>>> npp!1="nlon"
>>>>>>> ;;npp@lat2d=lat2d ; no ... only for graphics
>>>>>>> ;;npp@lon2d=lon2d ; no
>>>>>>> npp@coordinates="lat2d lon2d"
>>>>>>> npp@projection="sinusoidal"
>>>>>>
>>>>>> lat2d!0 = "nlat"
>>>>>> lat2d!1 = "nlon"
>>>>>> lon2d!0 = "nlat"
>>>>>> lon2d!1 = "nlon"
>>>>>>
>>>>>>>
>>>>>>> printVarSummary(npp)
>>>>>>>
>>>>>>> system("/bin/rm -f test.nc")
>>>>>>> out=addfile("test.nc","c")
>>>>>>> out->npp=npp
>>>>>>> out->lat2d=lat2d
>>>>>>> out->lon2d=lon2d
>>>>>>>
>>>>>>>
>>>>>>> On 8/30/11 12:05 PM, Michael Notaro wrote:
>>>>>>>> Thanks Dennis. Using hdfeof helps a lot and gives me the lats/lons now.
>>>>>>>>
>>>>>>>> Seems I am having another problem. After I read the hdfeos and retrieve
>>>>>>>> the lats/lons and data, I tried writing it to a file and always get segmentation fault.
>>>>>>>> What am I doing wrong?
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/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/csm/shea_util.ncl"
>>>>>>>> begin
>>>>>>>>
>>>>>>>> ; 1km annual npp modis on 1200x1200 tiles
>>>>>>>>
>>>>>>>> a=addfile("/Users/notaro/downloads/Y2000/MOD17A3.A2000366.h11v11.305.2009342143541.hdf.hdfeos","r")
>>>>>>>>
>>>>>>>> lat2d=a->GridLat_MOD_Grid_MOD17A3
>>>>>>>> lon2d=a->GridLon_MOD_Grid_MOD17A3
>>>>>>>>
>>>>>>>> npp=short2flt(a->Npp_1km_MOD_Grid_MOD17A3)
>>>>>>>>
>>>>>>>> npp!0="lat2d"
>>>>>>>> npp!1="lon2d"
>>>>>>>> npp@lat2d=lat2d
>>>>>>>> npp@lon2d=lon2d
>>>>>>>> npp@coordinates="lat2d lon2d"
>>>>>>>> npp@projection="sinusoidal"
>>>>>>>>
>>>>>>>> printVarSummary(npp)
>>>>>>>>
>>>>>>>> out=addfile("test.nc","c")
>>>>>>>> out->npp=npp
>>>>>>>> out->lat2d=lat2d
>>>>>>>> out->lon2d=lon2d
>>>>>>>>
>>>>>>>> end
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> Variable: npp
>>>>>>>> Type: float
>>>>>>>> Total Size: 5760000 bytes
>>>>>>>> 1440000 values
>>>>>>>> Number of Dimensions: 2
>>>>>>>> Dimensions and sizes: [lat2d | 1200] x [lon2d | 1200]
>>>>>>>> Coordinates:
>>>>>>>> Number Of Attributes: 7
>>>>>>>> lon2d : <ARRAY of 1440000 elements>
>>>>>>>> lat2d : <ARRAY of 1440000 elements>
>>>>>>>> _FillValue : 1e+20
>>>>>>>> coordinates : lat2d lon2d
>>>>>>>> hdfeos_name : Npp_1km
>>>>>>>> projection : sinusoidal
>>>>>>>> unsigned : True
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Aug 30, 2011, at 12:44 PM, Dennis Shea wrote:
>>>>>>>>
>>>>>>>>> Mike
>>>>>>>>>
>>>>>>>>> [1]
>>>>>>>>> %> ncl_filedump MOD17A3.A2010365.h00v08.305.2011029145245.hdf
>>>>>>>>>
>>>>>>>>> yields
>>>>>>>>>
>>>>>>>>> Variable: f
>>>>>>>>> Type: file
>>>>>>>>> filename: MOD17A3.A2010365.h00v08.305.2011029145245
>>>>>>>>> path: MOD17A3.A2010365.h00v08.305.2011029145245.hdf
>>>>>>>>> file global attributes:
>>>>>>>>> HDFEOSVersion : HDFEOS_V2.16
>>>>>>>>> StructMetadata_0 : GROUP=SwathStructure
>>>>>>>>> END_GROUP=SwathStructure
>>>>>>>>> GROUP=GridStructure
>>>>>>>>> GROUP=GRID_1
>>>>>>>>> GridName="MOD_Grid_MOD17A3"
>>>>>>>>> XDim=1200
>>>>>>>>> YDim=1200
>>>>>>>>> UpperLeftPointMtrs=(-20015109.354000,1111950.519667)
>>>>>>>>> LowerRightMtrs=(-18903158.834333,0.000000)
>>>>>>>>> Projection=GCTP_SNSOID
>>>>>>>>> [SNIP]
>>>>>>>>>
>>>>>>>>> [2]
>>>>>>>>> If you append the '.hdfeos' suffix, you will get that lat/lon values.
>>>>>>>>>
>>>>>>>>> %> ncl_filedump MOD17A3.A2010365.h00v08.305.2011029145245.hdf.hdfeos
>>>>>>>>>
>>>>>>>>> Variable: f
>>>>>>>>> Type: file
>>>>>>>>> filename: MOD17A3.A2010365.h00v08.305.2011029145245.hdf
>>>>>>>>> path: MOD17A3.A2010365.h00v08.305.2011029145245.hdf
>>>>>>>>> file global attributes:
>>>>>>>>> dimensions:
>>>>>>>>> YDim_MOD_Grid_MOD17A3 = 1200
>>>>>>>>> XDim_MOD_Grid_MOD17A3 = 1200
>>>>>>>>> variables:
>>>>>>>>> double GridLat_MOD_Grid_MOD17A3 ( YDim_MOD_Grid_MOD17A3,
>>>>>>>>> XDim_MOD_Grid_MOD17A3 )
>>>>>>>>> projection : (null)
>>>>>>>>> corners : ( 9.995833332438675, 9.995833332438675,
>>>>>>>>> 0.004166666666293064, 0.004166666666293064 )
>>>>>>>>> long_name : latitude
>>>>>>>>> units : degrees_north
>>>>>>>>>
>>>>>>>>> double GridLon_MOD_Grid_MOD17A3 ( YDim_MOD_Grid_MOD17A3,
>>>>>>>>> XDim_MOD_Grid_MOD17A3 )
>>>>>>>>> projection : (null)
>>>>>>>>> corners : ( 177.2297839752788, -172.6245418649626,
>>>>>>>>> -170.00416710093, -179.9958337931229 )
>>>>>>>>> long_name : longitude
>>>>>>>>> units : degrees_east
>>>>>>>>> [snip]
>>>>>>>>>
>>>>>>>>> [3] Use the hdfeos suffix. The names are different than using 'just' hdf
>>>>>>>>> however NCL does provide the georeferenced data coordinates.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> [4] You might look at:
>>>>>>>>>
>>>>>>>>> http://www.ncl.ucar.edu/Applications/binning.shtml
>>>>>>>>>
>>>>>>>>> Good luck
>>>>>>>>> D
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On 8/30/11 11:23 AM, Michael Notaro wrote:
>>>>>>>>>> I am trying to work with MODIS HDF data of the land variable 1-km annual NPP (net primary productivity).
>>>>>>>>>>
>>>>>>>>>> ftp://ftp.ntsg.umt.edu/pub/MODIS/Mirror/MOD17_Science_2010/MOD17A3/
>>>>>>>>>>
>>>>>>>>>> Each file is on a separate tile. I put the info below on a sample file's
>>>>>>>>>> info on the variable Npp_1km and the 1200x1200 grid and bounding rectangle.
>>>>>>>>>>
>>>>>>>>>> It is my understanding that I should grab the rectangles that overlap my study region,
>>>>>>>>>> which is Wisconsin. I think 4 of them overlap it.
>>>>>>>>>>
>>>>>>>>>> I can read in the data with addfile (or use ncl_convert2nc). But I am not sure
>>>>>>>>>> how to retrieve the proper lats/lons to plot the data correctly. I used short2flt also
>>>>>>>>>> when reading it. I am guessing that the data grid must not be a normal rectangle but
>>>>>>>>>> stretched as the satellite moves, since the range of the lats is 10 degrees but the
>>>>>>>>>> range of the lons is 60 degrees.
>>>>>>>>>>
>>>>>>>>>> Any suggestions on how to plot this in NCL is appreciated.
>>>>>>>>>>
>>>>>>>>>> Michael
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> YDim_MOD_Grid_MOD17A3 = 1200 ;
>>>>>>>>>> XDim_MOD_Grid_MOD17A3 = 1200 ;
>>>>>>>>>>
>>>>>>>>>> short Npp_1km(YDim_MOD_Grid_MOD17A3, XDim_MOD_Grid_MOD17A3) ;
>>>>>>>>>> Npp_1km:hdf_name = "Npp_1km" ;
>>>>>>>>>> Npp_1km:units = "kg_C/m^2" ;
>>>>>>>>>> Npp_1km:long_name = "MOD17A3 NPP--MODIS Gridded 1KM Annual Net Primary Productivity (NPP)" ;
>>>>>>>>>> Npp_1km:_FillValue = -1s ;
>>>>>>>>>> Npp_1km:valid_range = 0s, -36s ;
>>>>>>>>>> Npp_1km:calibrated_nt = 23 ;
>>>>>>>>>> Npp_1km:add_offset_err = 0. ;
>>>>>>>>>> Npp_1km:add_offset = 0. ;
>>>>>>>>>> Npp_1km:scale_factor_err = 0. ;
>>>>>>>>>> Npp_1km:scale_factor = 0.0001 ;
>>>>>>>>>>
>>>>>>>>>> " OBJECT = NORTHBOUNDINGCOORDINATE\n",
>>>>>>>>>> " NUM_VAL = 1\n",
>>>>>>>>>> " VALUE = 69.9999999937138\n",
>>>>>>>>>> " END_OBJECT = NORTHBOUNDINGCOORDINATE\n",
>>>>>>>>>> "\n",
>>>>>>>>>> " OBJECT = SOUTHBOUNDINGCOORDINATE\n",
>>>>>>>>>> " NUM_VAL = 1\n",
>>>>>>>>>> " VALUE = 59.9999999946118\n",
>>>>>>>>>> " END_OBJECT = SOUTHBOUNDINGCOORDINATE\n",
>>>>>>>>>> "\n",
>>>>>>>>>> " OBJECT = EASTBOUNDINGCOORDINATE\n",
>>>>>>>>>> " NUM_VAL = 1\n",
>>>>>>>>>> " VALUE = -119.983333303015\n",
>>>>>>>>>> " END_OBJECT = EASTBOUNDINGCOORDINATE\n",
>>>>>>>>>> "\n",
>>>>>>>>>> " OBJECT = WESTBOUNDINGCOORDINATE\n",
>>>>>>>>>> " NUM_VAL = 1\n",
>>>>>>>>>> " VALUE = -179.999999938971\n",
>>>>>>>>>> " END_OBJECT = WESTBOUNDINGCOORDINATE\n",
>>>>>>>>>> "\n",
>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>>> 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
>>>>>>> _______________________________________________
>>>>>>> ncl-talk mailing list
>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>