| reprojection [message #75071] |
Tue, 15 February 2011 21:42 |
gaurijyoti29
Messages: 10
Registered: April 2009
|
Junior Member |
|
|
I have extracted part of thid code from this group and tried to use
for my hdf file. However, I am not able to write the ouput file in UTM
projection. I request you all to kindly suggest me how to do it. The
code is as follows;
tx=24 ; tile horizontal number
ty=5 ; tile vertical number
utm_zone=43 ; UTM zone of the output projection (please choose the one
you need)
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
;reading in the HDF file:
hdfid=hdf_sd_start('H:\MOD16A2_MONTHLY
\MOD16A2.A2000M01.h24v05.105.2010358053227.hdf')
list=hdf_sd_varlist(hdfid)
print,list
hdf_sd_varread,hdfid,'ET_1km',ET_1km
nx=1200
ny=1200
; 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))
; transform into UTM coordinates
mapStruct = map_proj_init('integerized sinusoidal', center_longitude =
0., $
false_northing = 0., false_easting = 0.,
is_zones=86400.0, $
sphere_radius = 6371007.181)
limit = map_proj_inverse(lons,lats, map_structure=mapStruct)
xy=map_proj_forward(lons,lats,map_structure=mapStruct)
xutm=dblarr(nx,ny)
yutm=dblarr(nx,ny)
xutm[*]=xy[0,*]
yutm[*]=xy[1,*]
print,yutm/100000.0
print,xutm/100000.0
Would appreciate your cooperation.
Jyothi
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
