On 17 Feb., 21:59, David Fanning <n...@dfanning.com> wrote:
> Roman writes:
>> 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.
>
> Humm. Well, you seem surprised by this result, but nowhere
> in your code do you set up an integerized sinusoidal projection.
> You appear to be using a cylindrical projection. I'm pretty
> sure that if the data is in one projection, and the overlays
> are in another, that the chances of good alignment are pretty
> poor, indeed. :-)
>
> IDL does have an integerized sinusoidal projection, but
> you will have to use MAP_PROJ_INIT to set it up, and you
> will have to pass the map structure created by MAP_PROJ_INIT
> to Map_Continents and Map_Grid, so they will be able to
> draw on the map. I've never used this projection, but it has
> a number of 131. I am not sure it will be available in your
> version of IDL. It may have been added *after* IDL 5.6.
>
> Your code will look something like this:
>
> window, xsize=800, ysize=800
> position = [0.1, 0.1, 0.9, 0.9]
> mapStruct = Map_Proj_Init(131, ...)
> TVScale, image, Position=position, /KEEP
> Plot, mapStruct.uv_box[[0,2]], mapStruct.uv_box[[1,3]], $
> Pos=position, /NODATA, /NOERASE
> Map_Continents, MAP_STRUCTURE=mapStruct
> Map_Grid, MAP_STRUCTURE=mapStruct
>
> Cheers,
>
> David
>
> --
> David Fanning, Ph.D.
> Coyote's Guide to IDL Programming (www.dfanning.com)
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
Thanks for all your comments ! Below you find the way I did it
succesfully in IDL: Cheers, Roman
....
....
; 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,*]
xy=map_proj_forward(lonmin,latmin,map_structure=mapStruct)
xmin=xy[0]
ymin=xy[1]
xy=map_proj_forward(lonmax,latmax,map_structure=mapStruct)
xmax=xy[0]
ymax=xy[1]
print,'UTM boundary coordinates for the chosen area:
',xmin,xmax,ymin,ymax
;; the actual triangulation and reprojection goes here
triangulate,xutm,yutm,tr
dx=out_resolution*1000.
dy=out_resolution*1000.
....
....
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=[nodata,maxdata],colrange=[0,254],/
nomap,thick=2,xthick=2,ythick=2
Map_Continents,/COUNTRIES, /COASTS,/HIRES,mlinethick=2
Map_Grid,/box_axes,/clip_text,charsize=1.5,charthick=2,/HORI ZON
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
