Hi there
After trying out for several days to map (hundreds of) MODIS MOD14A2
products in IDL, I hope someone can help me here ! I went through
older comments on that but I really don't know where is my mistake !
The dataset and the boundaries of the continents/coast lines dont fitt
to each other. In the upper right corner yes but than it gets worse to
the lower right corner !
Thank you in advance for any help - cheers Roman
--------------------
Datasets are HDF files, 1200 x 1200 rows/columns, resoultion 1 km in
integerized sinusoidal projection.
Here is part of my code (I work with IDL 5.6 on windows, I do not have
ENVI nor the MODIS reprojection tool):
...
tx=18 ; tile horizontal number
ty=4 ; tile vertical number
utm_zone=32 ; UTM zone of the output projection (please choose the
one you need)
lonmin=0.0
lonmax=15.0
latmin=40.0
latmax=50.0
in_resolution=1.0 ; resolution of the input datasete in km (SIN
projection)
out_resolution=1.0; output resolution in km (UTM projection)
wx = round(43200./in_resolution) ; world x size
wy = round(21600./in_resolution) ; world y size
.....after reading in the HDF file:
HDF_SD_GETINFO, sds_id, DIMS=dims
HDF_SD_ATTRINFO,sds_id,HDF_SD_ATTRFIND
(sds_id,"valid_range"),DATA=hdfrange
HDF_SD_ATTRINFO,sds_id,HDF_SD_ATTRFIND
(sds_id,"_FillValue"),DATA=hdffill
HDF_SD_GETDATA, sds_id, array
hdffill=hdffill[0]
nx=dims[0]
ny=dims[1]
; create geolocation for specific tiles (lat/lon for each grid point)
lats=replicate(1.d,nx)#((90.d)-double(ty)/(18.d)*(180.d)-(di ndgen(ny)+.
5d)/double(wy)*(180.d))
lons=((dindgen(nx)+tx*nx)/double(wx)*(360.d)-(180.d))#((1.d) /(sin
((dindgen(ny)+ty*ny)*!dpi/double(wy))>1E-8))
map=map_proj_init(101,datum=8,zone=utm_zone)
xy=map_proj_forward(lons,lats,map_structure=map)
xutm=dblarr(nx,ny)
yutm=dblarr(nx,ny)
xutm[*]=xy[0,*]
yutm[*]=xy[1,*]
xy=map_proj_forward(lonmin,latmin,map_structure=map)
xmin=xy[0]
ymin=xy[1]
xy=map_proj_forward(lonmax,latmax,map_structure=map)
xmax=xy[0]
ymax=xy[1]
; the actual triangulation and reprojection goes here
triangulate,xutm,yutm,tr
dx=out_resolution*1000.
dy=out_resolution*1000.
dims=size(outdata,/dimensions)
ox=dims[0]
oy=dims[1]
; the output data [outdata], is in a regular UTM grid with grid bounds
[xmin, ymin, xmax, ymax] and should be plottet in a map
x1=latmax
x2=latmin
y1=lonmin
y2=lonmax
....create dummy file from the outdata file..^
sz=size(dummy) & nx=sz(1) & ny=sz(2) & mag=1
map_set,/grid,/continent,/HIRES, Position=[0.05, 0.05, 0.95, 0.80],$
charsize=1,LIMIT=[x2,y1,x1,y2]
imagemap,congrid(dummy,mag*nx,mag*ny), congrid(lats,mag*nx,mag*ny,/
interp), congrid(lons,mag*nx,mag*ny,/interp),$
PIXRANGE=[0,254],colrange=[0,254] ;Position=[0.05, 0.05, 0.95, 0.80],
PIXRANGE=[0,254],colrange=[0,254],thick=2,xthick=2,ythick=2
map_continents, /COASTS,/HIRES, mlinethick=2,/COUNTRIES
map_grid,/box_axes,/clip_text,charsize=1.5,charthick=2,/HORI ZON
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
