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