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