I have had bad luck trying to get map_set to do things right, or even
properly correlating the coordinate grid with the map features, though
I suppose it is possible that problem has been fixed by now.
I wrote my own, that reads the map files that IDL defines, and draws
things in a very simple manner.
After draw_map is called, user coordinates are identical to lat,lon
(degrees north and east) coordinates, so if you want to plot something
on the map, just do
PLOTS,LATITUDE_ARRAY,LONGITUDE_ARRAY
--------------------------CUT HERE---------------------
; Contains two routines, by Mitchell R Grunes:
; Nice_Tic: Make axis tic intervals nice.
; Draw_Map: Draw a map.
;----------------------------------------------------------- -----------------
; Routine to make axis tic intervals nice.
;----------------------------------------------------------- -----------------
pro Nice_Tic,xmin,xmax,absmin,absmax, xticks, xtickv
; ---------------INPUT---------------
; xmin and xmax=min,max plot values.
; absmin,absmax=boundaries for ticks--
; e.g., for latitude, should use -90,90
; xticks=maximum desired # of
; tick intervals (will be revised).
; ---------------OUTPUT--------------
; xticks, xtickv = the PLOT parameters
; of the same name.
;By Mitchell Grunes.
dx=double(xmax-xmin)/xticks ;Approximate interval length.
dxlog=alog10(dx)
dxlog=long(dxlog+100)-100
scale=10.d0^dxlog
dx=dx/scale ;Scaled to interval [1,10].
if dx lt 1.0d0 then begin
dx=1.0d0
endif else if dx lt 1.5d0 then begin
dx=1.5d0
endif else if dx le 2.0d0 then begin
dx=2.0d0
endif else if dx le 2.5d0 then begin
dx=2.5d0
endif else if dx le 3.0d0 then begin
dx=3.0d0
endif else if dx lt 5.0d0 then begin
dx=5.0d0
endif else begin
dx=10.d0
endelse
xmins=xmin/scale
if xmins lt absmin/scale then xmins=absmin/scale
if xmins ge 0 then begin
xmins= long( xmins/dx+1)*dx
endif else begin
xmins=-long(-xmins/dx )*dx
endelse
xmaxs=xmax/scale
if xmaxs gt absmax/scale then xmaxs=absmax/scale
if xmaxs ge 0 then begin
xmaxs= long( xmaxs/dx )*dx
endif else begin
xmaxs=-long(-xmaxs/dx )*dx
endelse
xticks=fix((xmaxs-xmins)/dx)
if xticks gt 30 then dx=dx*2 ; Do to IDL limits
xtickv=(indgen(xticks+1)*dx+xmins)*scale
end
;----------------------------------------------------------- -----------------
; Routine to Draw a map.
;----------------------------------------------------------- -----------------
pro Draw_Map,latmin,latmax,lonmin,lonmax,title
;Draw a map. It will then be possible to
; plot over the map using
; PLOTS,LON,LAT,NOCLIP=0
; (LON and LAT are arrays in degrees).
;latmin,latmax=minimum,maximum latitude (>0=north)
;lonmin,lonmax=minimum,maximum longitude (>0=east)
;title=plot title
;By Mitchell Grunes.
;Simplified from userlib procedure map_set.pro,
; because I could not figure out a way to get
; map_set.pro to give a simple rectangular
; projection, and because it seemed to be drawing
; things in the wrong places.
latmin2=latmin
latmax2=latmax
lonmin2=lonmin
lonmax2=lonmax
;Adjust aspect, lat and lon boundaries so that
;the lat and lon scales will be the same at map
;center. Approximately valid for default postscript
;in landscape mode.
aspect=8.3125/6.125/cos((latmin+latmax)/2.*!pi/180)
if (lonmax-lonmin) lt (latmax-latmin)*(aspect*.98) then begin
d=(latmax-latmin)*aspect-(lonmax-lonmin)
lonmin2=lonmin2-d/2
lonmax2=lonmax2+d/2
endif else if (latmax-latmin) lt (lonmax-lonmin)/(aspect*.98) then begin
d=(lonmax-lonmin)/aspect-(latmax-latmin)
latmin2=latmin2-d/2
latmax2=latmax2+d/2
endif
lonmin3=lonmin2 ;Keep within 360 degree range
lonmax3=lonmax2
if lonmax2-lonmin2 gt 360 then begin
avg=(lonmin2+lonmax3)/2.d0
lonmin3=avg-180.01d0
lonmax3=avg+180.01d0
endif
xticks=15
Nice_Tic,lonmin2,lonmax2,lonmin3,lonmax3, xticks, xtickv
yticks=20
Nice_Tic,latmin2,latmax2,-90, 90, yticks, ytickv
plot,[lonmin2,lonmin2,lonmax2,lonmax2,lonmin2], $
[latmin2,latmax2,latmax2,latmin2,latmin2],xstyle=1,ystyle=1, title=title, $
ticklen=1,xticks=xticks,yticks=yticks,xtickv=xtickv,ytickv=y tickv
close,1
openr,1,FILEPATH('supmap.dat',subdir = "maps"),/xdr,/stream
fbyte = [ 0, 71612L, 165096L]
nsegs = [ 283, 325, 594 ]
ij=2 ;0=course resolution map,1=U.S. only,2=All
point_lun, 1, fbyte(ij)
for i=1,nsegs(ij) do begin
npts = 0L
maxlat=0. & minlat=0. & maxlon=0. & minlon=0.
readu,1,npts,maxlat,minlat,maxlon,minlon
npts = npts / 2 ;# of points
xy = fltarr(2,npts)
readu,1,xy
lat = xy(0,*)
lon = xy(1,*)
; if stuff out of range, skip
; map segment.
if maxlat lt latmin2 or minlat gt latmax2 then goto,skip
if maxlon le lonmin3 then begin
minlon=minlon+360
maxlon=maxlon+360
lon=lon+360
endif
if minlon ge lonmax3 then begin
minlon=minlon-360
maxlon=maxlon-360
lon=lon-360
endif
if (maxlon lt lonmin3) or (minlon gt lonmax3) then goto,skip
plots, lon,lat,NOCLIP=0,color=3*!d.n_colors/9
empty
skip:
endfor
end
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
