David Fanning wrote:
> =
> I have to admit that I am confused about recent postings
> from St=E9phane Erard and Martin Schultz, who complain about
> having to fiddle with the Position parameters to get a
> Map_Set command to fit the proper Limits around an image.
>
> [...]
> =
> What I am unsure about is if the grid is going *exactly*
> where it is suppose to go. How would I know? You folks
> presumable have some kind of ground point or something
> that you know the exact latitude and longitude coordinate
> of. If I do the same thing with a data set I am more
> familiar with (say, the worldelv.dat data set in IDL), then
> it still appears to do the right thing.
> =
> In any case, I am all ears. :-)
> =
Ok, I'll try to make this problem a little clearer. I answer in the
Newsgroup since this issue probably interests other people.
Included is a small piece of source code that plots a map. The map is a
subset of the USGS shaded relief map of Mars, available for instance at
http://www-pdsimage.jpl.nasa.gov/PDS/. The subset is saved as
ptcarte.idl (it can be found at ftp://ftp.ias.fr/ias/erard), but any
image of this size will do to figure out the problem. You can see a gif
version of the (almost) final figure at
http://www.ias.fr/cdp/ISM/specomp.html
If you run the following code in two versions, you'll see a difference:
- as is, makes a correction to the map dimension, and the output is just
fine: the final product has the same number of pixel that the original
map file and the image fills the coordinates grid.
- If now you comment out the 6 lines beginning with del=3D... (ie use
map_set to plot the map without precautions) the output is different (at
least in versions 4 and 3.6). The frame now has the dimension of the
image (720 pixels in this case), so the image fits in the frame rather
than in the lat/long grid, which is 2% smaller.
The (spurious) transform is made in Map_uv_bound, a function in the
map_set.pro file (look for "fudge" in the text).
I understand that map_uv_bounds tries to plot a slightly larger frame so
we can see the axes clearly in the frame. Practically, it enlarges the
lat/long limits by 2% and plots these limits in the area reserved by the
keyword position in map_set. By doing so, it establishes a
correspondance between the area selected and lat/long limits that are 2%
in excess. In the exemple code, it sets the frame's right border to
long=3D 20+(20-(-70))*0.01.
If you don't set the map limits everything is fine - this explains why
there is no problem with the whole Earth map.
Now this problem is important, because:
- if you work with geographical or celestial images you know the
lat/long limits precisely, but you get a shrinked lat/long grid; this
precludes any precise position measurement of features of interest;
- if you try and superimpose another imaging data set on the first
image, IDL will use the coordinate system set by plot in map_uv_bound,
therefore corresponding to an area smaller than the first image. This
precludes precise superposition of geographical data sets, namely
construction of data cubes through the graphic interface. This is a
_major_ problem for a plotting/computing application (of course, this
will happen only if you plot the second data set with oplot or
polyfill).
My original question was about the reason why somebody took the pain to
make this change (modication made in march 1993, that should be version
3.6). Perhaps David is right in suggesting I'm optimistic in supposing
there is a reason, but still... I'm planning to adapt map_set to
planetary mapping, so I need to correct also problems like this one. If
it was a bug correction, I wouldn't like to resurect an old ghost.
Anybody has an idea?
Pro bid, d, PS=3DPS
; =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D =3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D =3D=3D=3D=3D=3D=
=3D=3D=3D
; example version from Mapirs5.pro. Only the map is plotted, no
additional data
; =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D =3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D =3D=3D=3D=3D=3D=
=3D=3D=3D
common colors,r_o,B_O,G_O,R,G,B
x=3Dfltarr(4)
y=3Dfltarr(4)
tab=3Dbytarr(1440,723)
lat=3Dfltarr(4)
lon=3Dfltarr(4)
!P.font=3D0
if keyword_set(ps) then begin
set_plot,'ps'
!P.font=3D0
coefPS=3Dmin([!D.X_vsize/850.,!D.Y_vsize/420.]) ;image on the whole shee=
t
device,filename=3D'Mapirs5.ps',/color,bits=3D8,/landscape,/B old,scale_fac=
tor=3D1.
endif else begin
window, /free, xsize=3D850, ysize=3D420
coefPS=3D1.
endelse
loadct, 0
l=3Dreplicate(0.7,16) ; define the rest of the color map for other data
s=3Dreplicate(1.,16)
h=3Dfindgen(16)*290./16.
h=3Dreverse(h)
h(0)=3D0.
L(0)=3D0.5
S(0)=3D1.
tvlct, H,L,S, 212,/HLS
tvlct, ro, gr,bl, /get
; following lines select the area of interest from the complete
map.
;e=3Da(440:799,280:439) ; lat: -20 to 20, long: -70 to 20 on the equator=
;a=3D1 ; (terrestrial convention for longitude). =
;d=3Drebin(e,720,320)
;save, filename=3D'ptcarte.idl', d,/XDR
; Reads the PDS map, enlarged, of the selected area
; (complete map is 1440 x 720). Sinusoidal projection
restore, 'ptcarte.idl' =
nbcol=3D211. ; first 201 values for the map
d=3Dbytscl(d,top=3Dnbcol)
sz=3Dsize(d)
px0 =3D20. ;shift a little
py0 =3D20.
print, sz, coefPS, !D.x_vsize,!D.y_vsize
qx=3D[px0,px0+sz(1)+0.5]/!D.x_vsize*coefPS ; Back to normal coord for
Map_set
qy=3D[py0,py0+sz(2)+0.5]/!D.y_vsize*coefPS
print, qx, qy
; these lines must be commented out to see what map_set does by
itself
del=3D(qx(1)-qx(0))*0.01 ;Compensate limits transform in
qx(0)=3Dqx(0)-del ; MAP_SET by a similar enlargement
qx(1)=3Dqx(1)+del ; of the plot.
del=3D(qy(1)-qy(0))*0.01
qy(0)=3Dqy(0)-del
qy(1)=3Dqy(1)+del
print, qx, qy, del
la=3D[-20,20]
lo=3D[-70, 20] ;At the equator
lax=3D(la(0)+la(1))/2.
lo1=3Dlo/cos(lax*!dtor) ;limits in long on top (latmax)
print, lo1, la
c0=3D[lax,lo1(0)] ; coord on the border centre =
c1=3D[la(1),(lo(0)+lo(1))/2.]
c2=3D[lax,lo1(1)]
c3=3D[la(0),(lo(0)+lo(1))/2.]
;print, lo1
tv, d, px0*coefPS,py0*coefPS,xsize=3Dsz(1)*coefPS, ysize=3Dsz(2)*coefPS
map_set, /sinu,/grid, limit=3D[c0,c1,c2,c3],/noerase, /label,title=3D$
'IRS-ISM overlaping observations of Mars!C',latlab=3D-72,lonlab=3D-22,$
londel=3D10.,latdel=3D10.,position=3D[qx(0),qy(0),qx(1),qy(1 )],$
glinethick=3D2.,glinestyle=3D0
maxi=3D0.5
mini=3D0.
;superpose another data set, no matter.
xyouts, 0,0, '!4'
colors=3Dreverse(findgen(15)+213) ; Couleurs albedo
citems=3Dstrarr(15)
citems(0)=3Dstring(format=3D'(F5.2)',maxi/32537./2.)
citems(7)=3Dstring(format=3D'(F5.2)',(maxi+mini)/32537./4.)
citems(14)=3Dstring(format=3D'(F5.2)',mini/32537./2.)
; use a routine from ASTRON
legend, citems,/fill,psym=3D8+intarr(15),colors=3Dcolors,char=3D1.,$
pos=3D[qx(1)+0.01,qy(1)],/normal
xyouts, qx(1)-0.015,qy(1)+0.02,/normal, 'ISM reflectance',charsize=3D0.8
!P.font=3D-1
if keyword_set(ps) then begin
device,/close
set_plot, 'x'
endif =
end
St=E9phane Erard
____________________________________________
Institut d'Astrophysique Spatiale
Universit=E9 Paris-Sud, bat. 121 F-91405 Orsay, France
tel : 33 1 69 85 86 41 fax : 33 1 69 85 86 75
e-mail : Erard@ias.fr
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
