1; MAP_PROJ_INIT 2; Documentation in progress 3; equivalent to the original procedure, except that instead of 40 or 4; so projections here we give access to all the PROJ projections (120+) 5; Note: the structure of the resulting mapstruct is not compatible with 6; IDL, i.e., one cannot use a GDL-defined mapstruct in IDL (if passed 7; within a save file), (but one can use an IDL-defined mapstruct in GDL). 8; March 2021: changed 'proj4...' to 'proj...' everywhere in GDL, except here where the structure tag is still 'proj4name' 9;rotates vector p1 around vector a by angle theta (degrees) 10function map_rotate3d, p1, a, theta 11 compile_opt idl2, hidden 12 ON_ERROR, 2 ; return to caller 13 deg2rad=!dpi/180d 14 st=sin(theta*deg2rad) 15 ct=cos(theta*deg2rad) 16 ; quaternion-derived rotation matrix 17 matrix=[[ a[0] * a[0] * (1 - ct) + ct, a[0] * a[1] * (1 - ct) - a[2]*st, a[0] * a[2] * (1 - ct) + a[1] * st],$ 18 [ a[1] * a[0] * (1 - ct) + a[2]*st, a[1] * a[1] * (1 - ct) + ct, a[1] * a[2] * (1 - ct) - a[0] * st],$ 19 [ a[2] * a[0] * (1 - ct) - a[1]*st, a[2] * a[1] * (1 - ct) + a[0]*st, a[2] * a[2] * (1 - ct) + ct ]] 20 ; multiply matrix vector 21 vector=[p1[0], p1[1], p1[2]] 22 rotated = matrix # vector 23 return, rotated 24 end 25 26pro map_proj_set_split, myMap 27 compile_opt idl2, hidden 28 ON_ERROR, 2 ; return to caller 29 ; v0 and u0 are the center lon & lat 30 sinlat = sin(myMap.v0) 31 coslat = cos(myMap.v0) 32 sinlon = sin(myMap.u0) 33 coslon = cos(myMap.u0) 34 xyzProjCenter = [coslon*coslat, sinlon*coslat, sinlat] 35; pole:lon & lat in rad, sin and cos of polar latitude, x,y,z 36; coordinates of pole. 37; SPLIT General case: if projection center is at (lon,lat) of 3D coords x,y,z and pole at xp,yp,zp 38; we have SPLIT=[lon, lat, CROSSP([x, y, z], [xp, yp, zp]), 0] 39 pole = myMap.pole[4:6] ;Location of pole 40 plane = CROSSP(xyzProjCenter, pole) 41 split=[myMap.p0lon, myMap.p0lat, plane, 0d] 42 MAP_CLIP_SET, MAP=myMap, SPLIT=split 43end 44 45; Ancillary procedure, compute limits (uv_box and ll_box) for full maps. Used only in initialization phase 46; of the projDefinitions.sav file 47pro gdl_compute_map_limits, myMap 48 49 compile_opt idl2, hidden 50 51 xmin=-1d & ymin=-1d & xmax=1d & ymax=1d 52 53 lonmin=-180.0d & lonmax=180.0d & latmin=-90.0d & latmax=90.0d 54 lonrange = lonmax - lonmin 55 latrange = latmax - latmin 56; is there another way (PROJ) to get ranges except brute force on a grid of possible points? 57; epsilon useful as projections are not precise (see #define EPS in PROJ c files: apparently < 1e-6) 58 epsx = 1d-6*ABS(lonrange) 59 epsy = 1d-6*ABS(latrange) 60 61 n_lons = 720 62 lons = [ lonmin + epsx, DINDGEN(n_lons)*(lonrange/(n_lons-1d)) + lonmin, lonmax - epsx] & n_lons += 2 63 64 n_lats = 360 65 lats = [latmin + epsy, DINDGEN(n_lats)*(latrange/(n_lats-1d)) + latmin, latmax - epsy, -epsy, epsy] & n_lats += 4 66 67 lons = reform(rebin(lons, n_lons, n_lats), 1, n_lons*n_lats) 68 lats = reform(rebin(transpose(lats), n_lons, n_lats), 1, n_lons*n_lats) 69 tmp = [temporary(lons), temporary(lats)] 70 71 xy = map_proj_forward(tmp, MAP=myMap) 72 ; Default if no points are valid is just the map_limits. 73 good = WHERE(FINITE(xy[0,*]) and FINITE(xy[1,*]), ngood) 74 75 if (ngood gt 0) then begin 76 xy = xy[*, good] 77 lonlat = tmp[*, good] 78 ; further check: are backprojected good points really close to original point? 79 tmp = MAP_PROJ_INVERSE(xy, MAP=myMap) 80 bad = WHERE(~FINITE(tmp[0,*]) or ~FINITE(tmp[1,*]), nbad) 81 if (nbad gt 0) then tmp[*, bad] = -9999 82 diff = ABS(lonlat - tmp) 83 diff[0,*] = diff[0,*] mod 360 ; Ignore 360 degre differences for longitude. 84 w = where(diff[0,*] gt 359.99, count) & if count gt 0 then diff[0,w] = 0 85 w = where((abs(tmp[0,*]) le 720) and (abs(tmp[1,*]) le 90) and (total(diff,1) lt 1d), count) 86 if (count gt 0) then begin ; Only those good during forward and inverse projection. 87 lonlat = lonlat[*, w] 88 xy = xy[*,w] 89 endif 90 xmin = min(xy[0,*], max=xmax) 91 ymin = min(xy[1,*], max=ymax) 92 lonmin = min(lonlat[0,*],poslonmin, max=lonmax, subscript_max=poslonmax) 93 latmin = min(lonlat[1,*],poslatmin, max=latmax, subscript_max=poslatmax) 94 endif 95 96 myMap.ll_box = [latmin, lonmin, latmax, lonmax] 97 myMap.uv_box = [xmin, ymin, xmax, ymax] 98end 99 100 101function map_proj_init, pindex, p4number=p4number, relaxed=relaxed, rotation=rotation, gctp=gctp, radians=radians, limit=passed_limit, ellipsoid=ellipsoid, semimajor_axis=semimajor_axis, semiminor_axis=semiminor_axis, sphere_radius=sphere_radius, center_azimuth=center_azimuth, datum=datum, clip=clip, gdl_precise=gdl_precise, check_proj4=chkprj4, _extra=extra 102 compile_opt idl2, hidden 103 ON_ERROR, 2 ; return to caller 104 105; NOTE: We are always "relaxed". 106; p4num bool indicates pindex is a number and refers to the internal PROJ table of PROJ properties line and not an IDL number for which an equivalent must be found 107 108; define limit 109if n_elements(passed_limit) lt 4 then limit=dblarr(4) else limit=double(passed_limit) 110 111; the common contains all relevant values after initialisation 112@gdlcommon_mapprojections_common 113deg2rad=!dpi/180d 114 115nkeys=n_tags(required) 116 117required_kw4=" +"+strlowcase(tag_names(required))+"=" 118 119nproj=n_elements(proj) 120 121sindex=pindex 122; find projection index, by index: 123if (N_ELEMENTS(pindex) le 0) then begin 124 index=where(proj.proj4name eq 'stere') ; stereo is default 125; this is a drawback: projection number is always an IDL number 126endif else if (SIZE(pindex, /TYPE) ne 7) then begin 127 if keyword_set(p4number) then sindex=proj[pindex].fullname else sindex=idl_ids[pindex] 128endif 129; now by name, as pindex is a string 130 131; grep best candidate name 132shortname = strupcase(strcompress(sindex, /REMOVE_ALL)) 133 ;defaults to IDL projection_names, check: 134w = strcmp(idl_ids,shortname,strlen(shortname)) & count=total(w) 135if count gt 1 then message, /noname, 'Ambiguous Projection abbreviation: ' + sindex 136if count eq 1 then begin 137 name4=idl_equiv[(where(w eq 1))[0]] 138 index=where(proj.proj4name eq name4, count) & if count eq 0 then message, 'Projection ' + sindex + ' apparently does not exist in Proj4 library, fixme.' 139endif else begin 140 ; next, PROJ projections 141 w = strcmp(compressed_ids1,shortname,strlen(shortname)) & count=total(w) 142 if count eq 0 then begin w = strcmp(compressed_ids2,shortname,strlen(shortname)) & count=total(w) & end ; alternative 143 if count eq 0 then message, /noname, 'Invalid Projection name: ' + sindex 144 if count eq 1 then index=where(w eq 1) else begin ; gt 1: ambiguous 145 message, /noname,/informational,'Ambiguous Projection name: ' + sindex 146 pos=where(w eq 1) 147 choices=compressed_ids1[pos] 148 lengths=strlen(choices) 149 match=min(lengths,j) 150 adopted=choices[j] 151 w = strcmp(compressed_ids1,adopted,strlen(compressed_ids1)) 152 pos=where(w eq 1) 153 index=pos[0] 154 endelse 155endelse 156 157; no more than nproj-1 158if index ge nproj then message, /noname, 'Invalid Projection number: ' + strtrim(sindex) 159 160; useful strings: 161; get base projection name! 162p4n=proj[index].proj4name 163; need to keep ony the real projection (PROJ) name if perchance there was additional commands already set in the name 164p4n=(strsplit(strtrim(p4n,2),' ',/extract))[0] 165 166; required parameters, filled. 167filled_required_parameter_string="" 168list_of_needed_params="" 169; optional parameters as passed 170filled_optional_parameter_string="" 171 172; create list of REQUIRED parameters 173required_opt="" 174for i=0,nkeys-1 do if required[index].(i) then required_opt+=required_kw4[i] 175if strlen(required_opt) gt 0 then begin 176 list_of_needed_params=strsplit(required_opt,' ',/extract) 177 n_required=n_elements(list_of_needed_params) 178endif else n_required=0 179 180; all optional elements in a string; add false northing easting that are 181; common for all projections. 182optional_opt=optional[index]+" x_0= y_0=" 183list_of_optional_params="+"+strsplit(optional_opt," ",/extract) 184n_optional = n_elements(list_of_optional_params) 185;print,list_of_optional_params 186 187;define type(s) of current projection 188property=proj_properties[index] 189if (property.EXIST eq 0) then message,'Unfortunately, projection '+shortname+ ' is flagged as absent. Please check MAP_INSTALL in GDL documentation.' 190conic=(property.CONIC eq 1) 191elliptic=(property.ELL eq 1) 192spheric=(property.SPH eq 1) 193cylindric=(property.CYL eq 1) 194azimuthal=(property.AZI eq 1) 195interrupted=(property.INTER eq 1) 196noRot=(property.NOROT eq 1) ; this indicates that, unexpectedly, a GOR cannot be tempted on this particular projection. 197; enable non-zero latitudes by use of general oblique rotation (GOR). 198; this is general but GOR adds uncertainties (and complexity=execution 199; time, especially near the poles). Besides, some projections crash 200; GDL when GOR is tempted. Finally, some projections do not need GOR 201; (they accept and use +lat_0=). Interrupted must not use GOR of course. 202; it can be forbidden by definition of NOROT 203rotPossible= ( ~noROT and ~azimuthal and ~conic and ~interrupted ) 204 205; GOR should not be tempted on projections that use +alpha and/or +lonc, as these 206; should be passed as arguments. 207if (n_required gt 0 or n_optional gt 0) then begin 208 w=WHERE(STRMATCH([list_of_needed_params,list_of_optional_params], '+alpha=') EQ 1, count) 209 if (count gt 0) then rotPossible=0 210endif 211 212south=0 213 214passed_params="" 215 216n_passed=0 217 218; enable abbreviated parameters 219if n_elements(extra) gt 0 then begin 220 221 passed_params = TAG_NAMES(extra) 222 223 ; instead of comparing passed_params directly to hash key, check possible abbreviations: 224 possible_params=tag_names(dictionary.tostruct()) 225 for i=0,n_tags(extra)-1 do begin 226 w = strcmp(possible_params,passed_params[i],strlen(passed_params[i])) & count=total(w) & j=(where(w eq 1))[0] 227 if count eq 1 then begin 228 passed_params[i]=possible_params[j] ; make passed_params normalized. 229 projkw=dictionary[possible_params[j]] 230 list_of_passed_params= (n_elements(list_of_passed_params) eq 0)? projkw : [list_of_passed_params,projkw] 231 passed_values= (n_elements(passed_values) eq 0)? strtrim(extra.(i),2) : [passed_values,strtrim(extra.(i),2)] 232 endif else if count gt 1 then message,"Ambiguous keyword abbreviation: "+passed_params[i] 233 endfor 234 n_passed=n_elements(list_of_passed_params) ; the only useful, recognized, ones, along with their passed_values. 235endif 236 237if n_passed gt 0 and n_passed ge n_required then begin 238; do we have required values for projection? 239 if n_required gt 0 then begin 240; needed params: get individual PROJ keywords, find if equivalent is 241; existing in passed_params. Based on equivalence list above. 242 243; sort in alphabetic order 244 s_needed=sort(list_of_needed_params) 245 ;create index table 246 tindex=(intarr(n_passed))-1; //-1 to further check 247 ; find all matches 248 for i=0,n_required-1 do begin 249 w=where(list_of_passed_params eq list_of_needed_params[i], count) & if count gt 0 then tindex[i]=w[0] ; we do not check duplicated entries, we take first. 250 endfor 251 ; tindex positive values need equal n_required here if all match were made 252 w=where(tindex ge 0, count, comp=absent) & if count ne n_required then begin 253 kwlist="" 254 ; reverse match: now we want needed_params that are not in passed_params. 255 for i=0,n_required-1 do begin 256 w=where(list_of_needed_params[i] eq list_of_passed_params, count) & if count eq 0 then kwlist+=yranoitcid[list_of_needed_params[i]]+" " 257 endfor 258 message,"Missing required parameters: "+kwlist 259 endif 260 ; populate required parameter list 261 for i=0,n_required-1 do begin 262 ; filter negative values for zone and set south 263 if list_of_needed_params[i] eq "+zone=" then begin 264 the_zone=fix(passed_values[tindex[i]]) 265 if the_zone lt 0 then begin 266 passed_values[tindex[i]]=strtrim(-1*the_zone,2) 267 south=1 268 endif 269 endif 270 filled_required_parameter_string+=" "+list_of_needed_params[i]+passed_values[tindex[i]] 271 endfor 272 endif 273 274 ; do we have optional values for projection? 275 if n_optional gt 0 then begin 276 ; sort in alphabetic order 277 ;create index table 278 tindex=(intarr(n_optional))-1 ; //-1 to further check 279 ; find all matches 280 for i=0,n_optional-1 do begin 281 w=where(list_of_passed_params eq list_of_optional_params[i], count) & if count gt 0 then tindex[i]=w[0] 282 endfor 283 ; w is the short list of optional passed values 284 w=where(tindex ge 0, count) & if count gt 0 then begin 285 ; populate optional parameter list 286 tindex=tindex[w] 287 for i=0,count-1 do begin 288 ; filter +n value between 0 and 1 only for fouc_s: 289 if list_of_passed_params[tindex[i]] eq "+n=" and p4n eq "fouc_s" then begin 290 n_val=fix(passed_values[tindex[i]]) 291 if (n_val lt 0 or n_val gt 1) then message,"Invalid parameter value ("+strtrim(n_val,2)+") for Foucaut Sinusoidal" 292 endif 293 filled_optional_parameter_string+=" "+list_of_passed_params[tindex[i]]+passed_values[tindex[i]] ; strtrim(extra.(tindex[i]),2) 294 endfor 295 endif 296 endif 297endif else begin ; or not... 298 if strlen(required_opt) gt 0 then message, "Absent (PROJ) parameter(s): "+required_opt 299endelse 300 301; main string (will need special treatment for rotation etc.) 302projcommand="+proj="+proj[index].proj4name+" " 303 304projoptions=filled_required_parameter_string+filled_optional_parameter_string 305 306; ok, projoptions contains all relevant AND permitted parameters. Try to assemble 307; all these into valid elements of !map. Up to now we have only 308; translated from idl to PROJ. now is time to interpret things a bit. 309 310; define useful defaults values 311p0lon = 0d ; center longitude 312p0lat = 0d ; center latitude 313 if n_elements(rotation) le 0 then rotation=0d ; rotation 314p1=0 ; locally used standard parallels p1 and p2 315p2=0 316satheight=0 317alpha=90 ; 318lonc=0 ; see f.e. oblique mercator in proj.org 319 320if strlen(projoptions) gt 0 then begin 321; convert projoptions to hash 322 s=strsplit(strtrim(projoptions,2),"= ",/extract) 323 x=where(strpos(s,'+') eq 0, comp=y) 324 a=hash(s[x],s[y]) 325; if a contains "+lon_0" this is p0lon, etc. 326 if a.HasKey("+lon_0") then p0lon=(a["+lon_0"]*1d)[0] 327 if a.HasKey("+lat_0") and rotPossible then begin 328 p0lat=(a["+lat_0"]*1d)[0] 329 endif 330 if a.HasKey("+lat_1") then p1=(a["+lat_1"]*1d)[0] 331 if a.HasKey("+lat_2") then p2=(a["+lat_2"]*1d)[0] 332 if a.HasKey("+h") then satheight=a["+h"]*1d 333 if a.HasKey("+alpha") then alpha=a["+alpha"]*1d 334 if a.HasKey("+lonc") then lonc=a["+lonc"]*1d 335endif 336 337 338; adjust ranges 339map_adjlon,p0lon 340p0lat= p0lat > (-89.999) & p0lat=p0lat < 89.999 341p1=p1 > (-89.999) & p1=p1<89.999 342p2=p2 > (-89.999) & p2=p2<89.999 343;if (p2 lt p1) then begin & tmp=p2 & p2=p1 & p1=p2 & end 344 345if (rotPossible) then begin 346 search_string='+lat_0=' 347; try a general oblique transformation 348 if n_passed ne 0 then begin 349 w=where(list_of_passed_params eq search_string, count) 350 if count gt 0 then begin ; try general oblique 351 p0lat=extra.(w[0]) 352 if p0lat ne 0 then begin 353 if p0lat gt 89.9 then p0lat = 89.9 ;take some precautions as PROJ is not protected!!! 354 if p0lat lt -89.9 then p0lat = -89.9 ; 355 ; compute pole of transformed projection 356 projcommand="+proj=ob_tran +o_proj="+proj[index].proj4name 357 ; remove '+lat_0=xxx +lon_0=xxx' from projoptions 358 a=strsplit(projoptions,"\+lat_0=[0-9.]*",/regex,/extract) 359 a=strsplit(projoptions,"\+lon_0=[0-9.]*",/regex,/extract) 360 projoptions=strjoin(a) 361 endif else rotPossible=0B 362 endif 363 endif else rotPossible=0B 364endif 365 366; for conic projections, although lat_0 is not in the list of 367; authorized parameters, it works, so we add it, as it is very 368; important to center the projection. 369if noRot eq 1 then p0lat=0 else begin 370 if (~rotPossible and (conic or spheric)) then begin ; 371 w=where(list_of_passed_params eq '+lat_0=', count) 372 if count gt 0 then begin 373 val=passed_values[w[0]] 374 p0lat=double(val) 375 projoptions+=" +lat_0="+val 376 endif 377 endif 378endelse 379 380; create a 999 !map 381myMap={!map} 382myMap.projection=999 383mymap.p[15]=index ;!useful for map_proj_info and unused apparently. 384myMap.p0lon = p0lon 385myMap.p0lat = p0lat 386myMap.u0 = p0lon * deg2rad 387myMap.v0 = p0lat * deg2rad 388;myMap.a = semimajor ; ellipsoid --> need table of correspondences! 389;myMap.e2 = e2 390 391myMap.rotation = rotation ; map rotation 392myMap.cosr=cos(rotation*deg2rad) 393myMap.sinr=sin(rotation*deg2rad) 394; pole is at +90 on meridian p0lon, eventually rotated by 395; center_azimuth 396pole_lon=p0lon 397pole_lat=p0lat+90 398; pole sines and xyz 399psinlat = 1d 400pcoslat = 0d 401psinlon = 0d 402pcoslon = 1d 403xyzpole = [pcoslon*pcoslat, psinlon*pcoslat, psinlat] 404; projection center sines and xyz 405csinlat = sin(myMap.v0) 406ccoslat = cos(myMap.v0) 407csinlon = sin(myMap.u0) 408ccoslon = cos(myMap.u0) 409xyzProjCenter = [ccoslon*ccoslat, csinlon*ccoslat, csinlat] 410 411if pole_lat gt 90 then begin 412 pole_lon+=180 & if pole_lon gt 180 then pole_lon-=360.0 413 pole_lat=180-pole_lat 414 psinlat = sin(pole_lat*deg2rad) 415 pcoslat = cos(pole_lat*deg2rad) 416 psinlon = sin(pole_lon*deg2rad) 417 pcoslon = cos(pole_lon*deg2rad) 418 xyzpole = [pcoslon*pcoslat, psinlon*pcoslat, psinlat] 419endif 420; rotate pole by center_azimuth if needed 421if keyword_set(center_azimuth) and rotPossible then begin 422 xyz=map_rotate3d(xyzpole,xyzprojcenter,center_azimuth) 423 pole_lon = atan(xyz[1], xyz[0])*!const.rtod 424 pole_lat = atan(xyz[2], sqrt(xyz[0]^2 + xyz[1]^2))*!const.rtod 425 xyzpole=xyz 426 psinlat = sin(pole_lat*deg2rad) 427 pcoslat = cos(pole_lat*deg2rad) 428endif 429 430; pole:lon & lat in rad, sin and cos of polar latitude, x,y,z 431myMap.pole=[pole_lon*deg2rad,pole_lat*deg2rad,psinlat,pcoslat,xyzpole] ; need to define myMap.pole BEFORE calling MAP_PROJ_SET_SPLIT 432 433; now that pole is computed correctly, add pole position to 434; generalized oblique 435if rotPossible then begin 436 projcommand+=" +o_alpha=90 +o_lat_c="+strtrim(p0lat,2)+" +o_lon_c=180" 437 if (p0lat gt 0) then projcommand+=" +lon_0="+strtrim(p0lon,2) else projcommand+=" +lon_0="+strtrim(180+p0lon,2) 438endif 439 440 441MAP_CLIP_SET, MAP=myMap, /RESET ;Clear clipping pipeline. 442 443; 1) get !map useful values 444; radius or ell or.. 445myMap.a=6370997.0d ; default 446myMap.e2=1 447hasRadius=0 448ellipticalusagerequired=(p4n eq "utm" or p4n eq "ups") 449if n_elements(sphere_radius) gt 0 then begin 450 hasRadius=1 451 myMap.a=sphere_radius[0] 452endif else begin 453 radius=6370997.0d ; myMap.a already ok 454endelse 455 456hasEll=0 457 458if n_elements(ellipsoid) gt 0 or n_elements(datum) gt 0 then begin ; case where myMap.a will be false! 459 hasEll=1 460 if n_elements(datum) gt 0 and n_elements(ellipsoid) eq 0 then ellipsoid=datum ; in case both are present. 461 if ~(SIZE(ellipsoid, /TYPE) eq 7) then begin ; must give name (index in list) 462 if ellipsoid gt 25 or ellipsoid lt 0 then message,"Invalid value for keyword ELLIPSOID: "+strtrim(ellipsoid,2) 463 ellipsoid=ellipsoid_proj[ellipsoid] 464 endif else begin 465 w = strcmp(strupcase(ellipsoid_idl),strupcase(ellipsoid),strlen(ellipsoid)) & count=total(w) 466 if count gt 0 then begin 467 pos=where(w gt 0) & ellipsoid=ellipsoid_proj[pos[0]] 468 endif 469 endelse 470endif 471 472; defining +a and +b overrides R and Ell 473hasDefinedEll=0 474if n_elements(SEMIMAJOR_AXIS) gt 0 then begin 475 hasRadius=0 476 hasEll=0 477 if ~n_elements(SEMIMINOR_AXIS) gt 0 then Message,"Keywords SEMIMAJOR_AXIS and SEMIMINOR_AXIS must both be supplied." 478 hasDefinedEll=1 479 myMap.a=SEMIMAJOR_AXIS 480 f=1.-(semiminor_axis/semimajor_axis) 481 myMap.e2 = 2*f-f^2 482endif 483 484; 2) split, clip.. 485; treat all interrupted cases separately (probably need to change !Map 486; pipeline size to accomodate for many-faceted projections: not done 487; but easy! 488; non-azimuthal projections: split at 180 degrees from map center, 489; then: 490; conic: will cut out pole region at some lat, and stop somewhere on 491; the other side, usually not far from the lat_2 if exists. 492; cylindric: should cut 'cylinder ends' ---> not done properly for 493; transverse? 494; azim: should cut somewhere: gnomonic cannot show one hemisphere, 495; other can, but will be very distorted. 496 497if (interrupted or p4n eq 'bipc' or p4n eq "bertin1953" or p4n eq "qsc" ) then begin 498 case p4n of 499 "bertin1953": BEGIN 500 theta = deg2rad * 16.5 501 MAP_CLIP_SET, map=myMap, SPLIT=[16.5, 42, -sin(theta), $ 502 cos(theta), 0., 0.] 503 myMap.p0lon=16.5 504 myMap.p0lat=42 505 END 506 "qsc": BEGIN 507 MAP_CLIP_SET, MAP=myMap, CLIP_PLANE=[xyzProjCenter, -0.5] 508 END 509 "bipc": BEGIN 510 splits = [-20,-110] + p0lon 511 512 for i=0,n_elements(splits)-1 do begin 513 theta = deg2rad * splits[i] 514 MAP_CLIP_SET, map=myMap, SPLIT=[splits[i], 0, -sin(theta), cos(theta), 0., 0.] 515 endfor 516 myMap.up_flags=1000 ; redefine "epsilon" due to precision problems in PROJ!!!! 517 518 END 519 520 "igh": BEGIN 521 splits = [-180, -40, -100, -20, 80] + p0lon + 180.0d 522 523 for i=0,n_elements(splits)-1 do begin 524 theta = deg2rad * splits[i] 525 MAP_CLIP_SET, map=myMap, SPLIT=[splits[i], 0, -sin(theta), cos(theta), 0., 0.] 526 endfor 527 myMap.up_flags=1000000 ; redefine "epsilon" due to precision problems in pr 528 529 END 530 531 "rhealpix": BEGIN 532 splits = [-3*45, -2*45, -45, 0, 45, 2*45, 3*45] + 180d + p0lon 533 534 for i=0,n_elements(splits)-1 do begin 535 theta = deg2rad * splits[i] 536 MAP_CLIP_SET, map=myMap, SPLIT=[splits[i], 0, -sin(theta), cos(theta), 0., 0.] 537 endfor 538 map_clip_set, map=myMap, SPLIT=[0,90,0,0,1d,-2d/3d] 539 map_clip_set, map=myMap, SPLIT=[0,-90,0,0,-1d,-2d/3d] 540 myMap.up_flags=10000 ; redefine "epsilon" due to precision problems in PROJ!!!! 541 END 542 543 "healpix": BEGIN 544 splits = [-3*45, -2*45, -45, 0, 45, 2*45, 3*45] + 180d + p0lon 545 546 for i=0,n_elements(splits)-1 do begin 547 theta = deg2rad * splits[i] 548 MAP_CLIP_SET, map=myMap, SPLIT=[splits[i], 0, -sin(theta), cos(theta), 0., 0.] 549 endfor 550 myMap.up_flags=10000 ; redefine "epsilon" due to precision problems in PROJ!!!! 551 END 552 ELSE: print,"Interrupted projection "+p4n+" is not yet properly taken into account in map_proj_init, please FIXME!" 553 endcase 554endif else begin ; not interrupted 555 if not azimuthal then begin 556 MAP_PROJ_SET_SPLIT,myMap ; standard cut 557; conics: clip around the poles 558 if conic then begin 559 ; apparently clipping is done 10 degrees above or below equator for 560 ; opposite hemisphere and at +75 degrees on same hemisphere unless the 561 ; standard parallels are not on the same side of equator, giving a 75 562 ; degree clip on both sides. 563 test1= (p1 ge 0.0) ? 1 : -1 564 test2= (p2 ge 0.0) ? 1 : -1 565 if (test1 eq test2) then begin 566 map_clip_set, map=myMap, clip_plane=[0,0,test1,sin(deg2rad*10.)] 567 map_clip_set, map=myMap, clip_plane=[0,0,-1*test2,sin(deg2rad*75.0)] 568 myMap.p[13]=-1*test1*10. ; use it to store this value, see map_grid, map_horizon 569 myMap.p[14]=test2*75. ; use it to store this value, see map_grid, map_horizon 570 endif else begin 571 map_clip_set, map=myMap, clip_plane=[0,0,1,sin(deg2rad*75.0)] 572 map_clip_set, map=myMap, clip_plane=[0,0,-1,sin(deg2rad*75.0)] 573 myMap.p[13]=75. ; use it to store this value, see map_grid, map_horizon 574 myMap.p[14]=-75. ; use it to store this value, see map_grid, map_horizon 575 endelse 576 endif else if cylindric then begin ; mercator et: clip 10 degrees from poles. 577 ;alpha is the rotation angle for 578 ;oblique cylindric. alpha=0 for 579 ;transverse projections 580 val=80 581 ; remove poles 582 map_clip_set, map=myMap, clip_plane=[myMap.pole[4:6],sin(deg2rad*val)] 583 map_clip_set, map=myMap, clip_plane=[-1*myMap.pole[4:6],sin(deg2rad*val)] 584 myMap.p[13]=val ; use it to store this value, see map_grid, map_horizon 585 myMap.p[14]=-val ; use it to store this value, see map_grid, map_horizon 586 endif 587 endif else begin ; azim projs. 588 val=-1d-8 589 case p4n of 590 "tpers": if satheight eq 0 then val=-0.5d else val=-1.01d /(1+satheight/myMap.a) 591 "nsper": if satheight eq 0 then val=-0.5d else val=-1.01d /(1+satheight/myMap.a) 592 "gnom": val=-0.5d 593 else: val=-0.49999d 594 endcase 595 MAP_CLIP_SET, MAP=myMap, CLIP_PLANE=[xyzProjCenter, val] 596 myMap.p[14]=val ; use it to store this value, see map_grid, map_horizon 597 endelse ; end azim projs 598endelse ; not interrupted 599 600if ellipticalusagerequired then projOptions+=' +ellps=GRS80 ' 601if south then projOptions+=' +south' 602; finalize projection to be used in finding limits: 603myMap.up_name=projcommand+" "+projOptions 604 605if (hasDefinedEll) then begin 606myMap.up_name+=" +a="+strtrim(SEMIMAJOR_AXIS[0],2)+" +b="+strtrim(SEMIMINOR_AXIS[0],2) 607endif else if (hasRadius and ~ellipticalusagerequired) then begin 608myMap.up_name+=" +R="+strtrim(SPHERE_RADIUS[0],2) 609endif else if (hasEll) then begin 610myMap.up_name+=" +ell="+ellipsoid 611endif 612 613if keyword_set(chkprj4) then print,myMap.up_name 614; 3) Set LIMITs and clip. 615 616if keyword_set(gdl_precise) then gdl_compute_map_limits, myMap else gdl_set_map_limits, myMap, limit 617 618; 4) transform 619MAP_CLIP_SET, MAP=myMap, /transform ;apply transform 620 621return,myMap 622end 623 624; in case one wants to add new projections to the ../resource/maps/projections.csv file. 625pro map_proj_auxiliary_read_csv 626 compile_opt idl2 627 628 ON_ERROR, 2 ; return to caller 629 restore,"csv.sav" 630 nproj=n_elements(csv_proj.field1) 631 names={PROJ4NAME:"",FULLNAME:"",OTHERNAME:""} 632 proj_property={EXIST:1B,SPH:0B,CONIC:0B,AZI:0B,ELL:0B,CYL:0B,MISC:0B,NOINV:0B,NOROT:0B,INTER:0b} ; note uppercase 633; fill the sorted, uniq, list of REQUIRED values 634 t=strtrim(csv_proj.field5,2) & s=strjoin(t) & t=strsplit(s,"= ",/extract) 635 q=t[sort(t)] & required_template_list=q[uniq(q)] 636 for i=0,n_elements(required_template_list)-1 do map_struct_append, required_template, required_template_list[i], 0b 637 638 proj=replicate(names,nproj) 639 proj.PROJ4NAME=csv_proj.FIELD1 640 proj.FULLNAME=csv_proj.FIELD2 641 proj.OTHERNAME=csv_proj.FIELD3 642 643 proj_properties=replicate(proj_property,nproj) 644 csv_proj.field4=strupcase(strcompress(csv_proj.field4,/remove_all)) ; note uppercase 645 ntags=n_tags(proj_property) 646 tname=strupcase(tag_names(proj_property)) 647 for i=0,ntags-1 do proj_properties[ WHERE(STRMATCH(csv_proj.field4, '*'+tname[i]+'*') EQ 1)].(i)=1 648 649 required=replicate(required_template,nproj) 650 ntags=n_tags(required_template) 651 tname=strlowcase(tag_names(required_template)) ; ALL LOWCASE in table FOR THE MOMENT. WARNING if NOT!!! 652 ; to avoid problems, each element of csv_proj.field5, which is a 653 ; serie of strings like "plat_0= plon_0= phdg_0=" must be comparable 654 ; only with same tname (ex: plat_0) and not a subset tname such as 655 ; "lat_0". this implies to add at least one whitespace at the 656 ; beginning: 657 newfield5=" "+csv_proj.field5 658for i=0,ntags-1 do begin & w=WHERE(STRMATCH(newfield5, '* '+tname[i]+'=*', /FOLD_CASE) EQ 1, count) & if (count gt 0) then required[w].(i)=1 & end 659 660; optional 661 optional=csv_proj.field6 662; 663proj_limits=reform(replicate(-1.0d,4*nproj),4,nproj) 664proj_scale=dblarr(nproj) 665 666; save once to have map_proj_init work. 667save,filen="projDefinitions.sav",proj,proj_properties,required,optional,proj_scale,proj_limits 668 669; now compute default limits [-180..180, -90..90] for all projections. 670; the idea is to call all the projections with all 'possible' 671; parameters in order to have only unexisting projections that cause a 672; (trapped) error. As the projection has been set up, the uv box is 673; the one computed in map_proj_init using a brute force 674; method. Obvioulsy this could be made more exact if the uv_box was 675; part of the database, but imho this small amount of work will rebuke 676; everybody. An other option would be to use Proj functions to get all 677; the needed information, I've not looked into that. 678; call proj_init for uv_box approximate calculation... 679for i=0,nproj-1 do begin 680 catch,absent 681 if absent ne 0 then begin 682 print,'i was',i,' projection was ',proj[i].PROJ4NAME 683 proj_properties[i].exist=0b 684 continue 685 endif 686 myMap=map_proj_init(/gdl_precise, i,/p4num,alpha=0.0001,height=1,standard_parall=30,standard_par1=50,standard_par2=-45,sat_tilt=45,true_scale_latitude=12,lat_3=13,HOM_LONGITUDE1=1,HOM_LONGITUDE2=80,LON_3=120,OEA_SHAPEN=1, OEA_SHAPEM=1,SOM_LANDSAT_NUMBER=2, SOM_LANDSAT_PATH=22, ZONE=28) ; uses ellps=wgs84 by default. 687 proj_scale[i]=abs(myMap.uv_box[2]-myMap.uv_box[0]) ; number of ellipsoid meters in uv_box 688endfor 689; proj_scale, only on existing projections. 690for i=0,nproj-1 do begin 691 catch,absent 692 if absent ne 0 then begin 693 print,'(known?) problem with projection '+proj[i].PROJ4NAME 694 continue 695 endif 696 697 if proj_properties[i].exist eq 1 then begin 698 myMap=map_proj_init(/gdl_precise,i,/p4num,alpha=0.0001,sphere=1,height=1,standard_parall=30,standard_par1=50,standard_par2=-45,sat_tilt=45,true_scale_latitude=12,lat_3=13,HOM_LONGITUDE1=1,HOM_LONGITUDE2=80,LON_3=120,OEA_SHAPEN=1, OEA_SHAPEM=1,SOM_LANDSAT_NUMBER=2, SOM_LANDSAT_PATH=22, ZONE=28, /check) 699 proj_limits[*,i]=myMap.uv_box ; normalized uv_box 700 endif 701endfor 702 703 704save,filen="projDefinitions.sav",proj,proj_properties,required,optional,proj_scale,proj_limits 705end 706 707;soffice --headless --convert-to csv projections.ods 708;csv_proj=read_csv("projections.csv",n_table=1) & save,csv_proj,filename="csv.sav" & exit 709; gdl 710; .compile map_proj_init.pro 711; MAP_PROJ_AUXILIARY_READ_CSV 712; exit 713 714