1; ancillary MAP* procedure, set limits (uv_box and ll_box) for maps.
2pro gdl_set_map_limits, myMap, limits
3
4  compile_opt idl2, hidden
5
6; map_limits may have 4 or 8 elements.
7
8  if n_elements(limits) ne 8 and n_elements(limits) ne 4 then  message, 'Map limit must have 4 or 8 points'
9  map_limits=limits ; do not change limits
10  xmin=-1d & ymin=-1d & xmax=1d & ymax=1d
11
12  if n_elements(map_limits) eq 4 then begin
13; filter map_limits to [-180,+180] in both directions however keeping orientation
14     lonmin=map_limits[1] & lonmax=map_limits[3] & latmin=map_limits[0] & latmax=map_limits[2]
15; longitudes: beautify
16     map_adjlon,lonmin
17     map_adjlon,lonmax
18     if (lonmin gt lonmax) then  begin temp=lonmin & lonmin=lonmax & lonmax=temp & end
19     if (lonmin eq lonmax) then  begin lonmin=-180 & lonmax=180 & endif
20; latitudes: beautify
21     latmin = -90 > latmin < 90
22     latmax = -90 > latmax < 90
23     if (latmin gt latmax) then  begin temp=latmin & latmin=latmax & latmax=temp & end
24     if (latmin eq latmax) then  begin latmin=-90 & latmax=90 & endif
25
26     map_limits=[latmin, lonmin, latmax, lonmax]
27     if ( map_limits[0] eq  map_limits[2])  then message,/informational,"GDL_SET_MAP_LIMITS: Warning, MAP limits are invalid."
28; range of values
29       lonrange = lonmax - lonmin
30       latrange = latmax - latmin
31; is there another way (PROJ) to get ranges except brute force on a grid of possible points?
32; epsilon useful as projections are not precise (see #define EPS in PROJ c files: apparently < 1e-6)
33       epsx = 1d-6*ABS(lonrange)
34       epsy = 1d-6*ABS(latrange)
35
36       n_lons = 90
37       lons = [ lonmin + epsx, DINDGEN(n_lons)*(lonrange/(n_lons-1d)) + lonmin, lonmax - epsx] & n_lons += 2
38
39       n_lats = 45
40       lats =  [latmin + epsy, DINDGEN(n_lats)*(latrange/(n_lats-1d)) + latmin, latmax - epsy] & n_lats += 2
41     if ((latmin lt 0) && (latmax gt 0)) then begin lats = [lats, -epsy, epsy] &  n_lats += 2 & end
42
43       lons = reform(rebin(lons, n_lons, n_lats), 1, n_lons*n_lats)
44       lats = reform(rebin(transpose(lats), n_lons, n_lats), 1, n_lons*n_lats)
45       tmp = [temporary(lons), temporary(lats)]
46
47       xy = map_proj_forward(tmp, MAP=myMap)
48                                ; Default if no points are valid is just the map_limits.
49       good = WHERE(FINITE(xy[0,*]) and FINITE(xy[1,*]), ngood)
50
51       if (ngood gt 0) then begin
52         xy = xy[*, good]
53         lonlat = tmp[*, good]
54                                ; further check: are backprojected good points really close to original point?
55         tmp = MAP_PROJ_INVERSE(xy, MAP=myMap)
56         bad = WHERE(~FINITE(tmp[0,*]) or ~FINITE(tmp[1,*]), nbad)
57         if (nbad gt 0) then tmp[*, bad] = -9999
58         diff = ABS(lonlat - tmp)
59         diff[0,*] = diff[0,*] mod 360 ; Ignore 360 degre differences for longitude.
60         w = where(diff[0,*] gt 359.99, count) & if count gt 0 then diff[0,w] = 0
61         w = where((abs(tmp[0,*]) le 720) and (abs(tmp[1,*]) le 90) and (total(diff,1) lt 1d), count)
62         if (count gt 0) then begin ; Only those good during forward and inverse projection.
63           lonlat = lonlat[*, w]
64           xy = xy[*,w]
65         endif
66         xmin = min(xy[0,*], max=xmax)
67         ymin = min(xy[1,*], max=ymax)
68         lonmin = min(lonlat[0,*], max=lonmax)
69         latmin = min(lonlat[1,*], max=latmax)
70       endif
71
72  endif else begin
73; easy: convert and pray.
74; 8 point limit as in [latLeft,lonLeft, latTop, lonTop, LatRight, lonRight, LatBottom, LonBottom]
75;
76     lons=map_limits[[1,3,5,7]] & lats=map_limits[[0,2,4,6]]
77     xy=map_proj_forward(lons,lats,map=myMap)
78     good = WHERE(FINITE(xy[0,*]) and FINITE(xy[1,*]), ngood)
79     if (ngood ge 2 ) then begin
80        xmin = min(xy[0,*], max=xmax)
81        ymin = min(xy[1,*], max=ymax)
82        lonmin = min(lons, max=lonmax)
83        latmin = min(lats, max=latmax)
84     endif else message, 'Unmappable limit point(s) in LIMIT keyword'
85  endelse
86
87                                ; Fill in map structure.
88  myMap.ll_box = [latmin, lonmin, latmax, lonmax]
89  myMap.uv_box = [xmin, ymin, xmax, ymax]
90
91end
92
93