1 /******************************************************************************
2  *
3  * Purpose: Translation from ILWIS coordinate system information.
4  * Author:   Lichun Wang, lichun@itc.nl
5  *
6  ******************************************************************************
7  * Copyright (c) 2004, ITC
8  * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
9  *
10  * Permission is hereby granted, free of charge, to any person obtaining a
11  * copy of this software and associated documentation files (the "Software"),
12  * to deal in the Software without restriction, including without limitation
13  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  * and/or sell copies of the Software, and to permit persons to whom the
15  * Software is furnished to do so, subject to the following conditions:
16  *
17  * The above copyright notice and this permission notice shall be included
18  * in all copies or substantial portions of the Software.
19  *
20  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  * DEALINGS IN THE SOFTWARE.
27  ****************************************************************************/
28 #include "cpl_conv.h"
29 #include "ilwisdataset.h"
30 
31 #include <string>
32 
33 CPL_CVSID("$Id: ilwiscoordinatesystem.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
34 
35 namespace GDAL {
36 
37 typedef struct
38 {
39     const char  *pszIlwisDatum;
40     const char  *pszWKTDatum;
41     int   nEPSGCode;
42 } IlwisDatums;
43 
44 typedef struct
45 {
46     const char  *pszIlwisEllips;
47     int   nEPSGCode;
48     double semiMajor;
49     double invFlattening;
50 } IlwisEllips;
51 
52 static const IlwisDatums iwDatums[] =
53 {
54     { "Adindan", "Adindan", 4201 },
55     { "Afgooye", "Afgooye", 4205 },
56     // AGREF --- skipped
57     { "Ain el Abd 1970", "Ain_el_Abd_1970", 4204 },
58     { "American Samoa 1962", "American_Samoa_1962", 4169 },
59     //Anna 1 Astro 1965 --- skipped
60     { "Antigua Island Astro 1943", "Antigua_1943", 4601 },
61     { "Arc 1950", "Arc_1950", 4209 },    //Arc 1950
62     { "Arc 1960", "Arc_1960", 4210 },    //Arc 1960
63     //Ascension Island 1958
64     //Astro Beacon E 1945
65     //Astro DOS 71/4
66     //Astro Tern Island (FRIG) 1961
67     //Astronomical Station 1952
68     { "Australian Geodetic 1966", "Australian_Geodetic_Datum_1966", 4202 },
69     { "Australian Geodetic 1984", "Australian_Geodetic_Datum_1984", 4203 },
70     //Ayabelle Lighthouse
71     //Bellevue (IGN)
72     { "Bermuda 1957", "Bermuda_1957", 4216 },
73     { "Bissau", "Bissau", 4165 },
74     { "Bogota Observatory  (1975)", "Bogota", 4218 },
75     { "Bukit Rimpah", "Bukit_Rimpah", 4219 },
76     //Camp Area Astro
77     { "Campo Inchauspe", "Campo_Inchauspe", 4221 },
78     //Canton Astro 1966
79     { "Cape", "Cape", 4222 },
80     //Cape Canaveral
81     { "Carthage", "Carthage", 4223 },
82     { "CH1903", "CH1903", 4149 },
83     //Chatham Island Astro 1971
84     { "Chua Astro", "Chua", 4224 },
85     { "Corrego Alegre", "Corrego_Alegre", 4225 },
86     //Croatia
87     //D-PAF (Orbits)
88     { "Dabola", "Dabola_1981", 4155 },
89     //Deception Island
90     //Djakarta (Batavia)
91     //DOS 1968
92     //Easter Island 1967
93     //Estonia 1937
94     { "European 1950 (ED 50)", "European_Datum_1950", 4154 },
95     //European 1979 (ED 79
96     //Fort Thomas 1955
97     { "Gan 1970", "Gandajika_1970", 4233 },
98     //Geodetic Datum 1949
99     //Graciosa Base SW 1948
100     //Guam 1963
101     { "Gunung Segara", "Gunung_Segara", 4613 },
102     //GUX 1 Astro
103     { "Herat North", "Herat_North", 4255 },
104     //Hermannskogel
105     //Hjorsey 1955
106     //Hong Kong 1963
107     { "Hu-Tzu-Shan", "Hu_Tzu_Shan", 4236 },
108     //Indian (Bangladesh)
109     //Indian (India, Nepal)
110     //Indian (Pakistan)
111     { "Indian 1954", "Indian_1954", 4239 },
112     { "Indian 1960", "Indian_1960", 4131 },
113     { "Indian 1975", "Indian_1975", 4240 },
114     { "Indonesian 1974", "Indonesian_Datum_1974", 4238 },
115     //Ireland 1965
116     //ISTS 061 Astro 1968
117     //ISTS 073 Astro 1969
118     //Johnston Island 1961
119     { "Kandawala", "Kandawala", 4244 },
120     //Kerguelen Island 1949
121     { "Kertau 1948", "Kertau", 4245 },
122     //Kusaie Astro 1951
123     //L. C. 5 Astro 1961
124     { "Leigon", "Leigon", 4250 },
125     { "Liberia 1964", "Liberia_1964", 4251 },
126     { "Luzon", "Luzon_1911", 4253 },
127     //M'Poraloko
128     { "Mahe 1971", "Mahe_1971", 4256 },
129     { "Massawa", "Massawa", 4262 },
130     { "Merchich", "Merchich", 4261 },
131     { "MGI (Hermannskogel)", "Militar_Geographische_Institute",4312 },
132     //Midway Astro 1961
133     { "Minna", "Minna", 4263 },
134     { "Montserrat Island Astro 1958", "Montserrat_1958", 4604 },
135     { "Nahrwan", "Nahrwan_1967", 4270 },
136     { "Naparima BWI", "Naparima_1955", 4158 },
137     { "North American 1927 (NAD 27)", "North_American_Datum_1927", 4267 },
138     { "North American 1983 (NAD 83)", "North_American_Datum_1983", 4269 },
139     //North Sahara 1959
140     { "NTF (Nouvelle Triangulation de France)", "Nouvelle_Triangulation_Francaise", 4807 },
141     //Observatorio Meteorologico 1939
142     //Old Egyptian 1907
143     { "Old Hawaiian", "Old_Hawaiian", 4135 },
144     //Oman
145     //Ordnance Survey Great Britain 1936
146     //Pico de las Nieves
147     //Pitcairn Astro 1967
148     //Point 58
149     { "Pointe Noire 1948", "Pointe_Noire", 4282 },
150     { "Porto Santo 1936", "Porto_Santo",4615 },
151     //Potsdam (Rauenburg)
152     { "Potsdam (Rauenburg)", "Deutsches_Hauptdreiecksnetz", 4314 },
153     { "Provisional South American 1956", "Provisional_South_American_Datum_1956", 4248 },
154     //Provisional South Chilean 1963
155     { "Puerto Rico", "Puerto_Rico", 4139 },
156     { "Pulkovo 1942", "Pulkovo_1942", 4178 },
157     //{ "Qatar National", "Qatar_National_Datum_1995", 4614 },
158     { "Qornoq", "Qornoq", 4287 },
159     { "Puerto Rico", "Puerto_Rico", 4139 },
160     //Reunion
161     { "Rome 1940", "Monte_Mario", 4806 },
162     { "RT90", "Rikets_koordinatsystem_1990", 4124 },
163     { "Rijks Driehoeksmeting", "Amersfoort", 4289 },
164     { "S-42 (Pulkovo 1942)", "Pulkovo_1942", 4178 },
165     //{ "S-JTSK", "Jednotne_Trigonometricke_Site_Katastralni", 4156 },
166     //Santo (DOS) 1965
167     //Sao Braz
168     { "Sapper Hill 1943", "Sapper_Hill_1943", 4292 },
169     { "Schwarzeck", "Schwarzeck", 4293 },
170     { "Selvagem Grande 1938", "Selvagem_Grande", 4616 },
171     //vSGS 1985
172     //Sierra Leone 1960
173     { "South American 1969", "South_American_Datum_1969", 4291 },
174     //South Asia
175     { "Tananarive Observatory 1925", "Tananarive_1925", 4297 },
176     { "Timbalai 1948", "Timbalai_1948", 4298 },
177     { "Tokyo", "Tokyo", 4301 },
178     //Tristan Astro 1968
179     //Viti Levu 1916
180     { "Voirol 1874", "Voirol_1875", 4304 },
181     //Voirol 1960
182     //Wake Island Astro 1952
183     //Wake-Eniwetok 1960
184     { "WGS 1972", "WGS_1972", 4322 },
185     { "WGS 1984", "WGS_1984", 4326 },
186     { "Yacare", "Yacare", 4309 },
187     { "Zanderij", "Zanderij", 4311 },
188     { nullptr, nullptr, 0 }
189 };
190 
191 static const IlwisEllips iwEllips[] =
192 {
193     { "Sphere", 7035, 6371007, 0.0  },  //rad 6370997 m (normal sphere)
194     { "Airy 1830", 7031, 6377563.396, 299.3249646 },
195     { "Modified Airy", 7002, 6377340.189, 299.3249646 },
196     { "ATS77", 7204, 6378135.0, 298.257000006 },
197     { "Australian National", 7003, 6378160, 298.249997276 },
198     { "Bessel 1841", 7042, 6377397.155, 299.1528128},
199     { "Bessel 1841 (Japan By Law)", 7046 , 6377397.155, 299.152815351 },
200     { "Bessel 1841 (Namibia)", 7006, 6377483.865, 299.1528128 },
201     { "Clarke 1866", 7008, 6378206.4, 294.9786982 },
202     { "Clarke 1880", 7034, 6378249.145, 293.465 },
203     { "Clarke 1880 (IGN)", 7011, 6378249.2, 293.466 },
204     // FIXME: D-PAF (Orbits) --- skipped
205     // FIXME: Du Plessis Modified --- skipped
206     // FIXME: Du Plessis Reconstituted --- skipped
207     { "Everest (India 1830)", 7015, 6377276.345, 300.8017 },
208     // Everest (India 1956) --- skipped
209     // Everest (Malaysia 1969) --- skipped
210     { "Everest (E. Malaysia and Brunei)", 7016, 6377298.556, 300.8017 },
211     { "Everest (Malay. and Singapore 1948)", 7018, 6377304.063, 300.8017 },
212     { "Everest (Pakistan)", 7044, 6377309.613, 300.8017 },
213     // Everest (Sabah Sarawak) --- skipped
214     // Fischer 1960 --- skipped
215     // Fischer 1960 (Modified) --- skipped
216     // Fischer 1968 --- skipped
217     { "GRS 80", 7019, 6378137, 298.257222101  },
218     { "Helmert 1906", 7020, 6378200, 298.3 },
219     // Hough 1960 --- skipped
220     { "Indonesian 1974", 7021, 6378160, 298.247 },
221     { "International 1924", 7022, 6378388, 297 },
222     { "Krassovsky 1940", 7024, 6378245, 298.3 },
223     // New_International 1967
224     // SGS 85
225     // South American 1969
226     // WGS 60
227     // WGS 66
228     { "WGS 72", 7020, 6378135.0, 298.259998590  },
229     { "WGS 84", 7030, 6378137, 298.257223563 },
230     { nullptr, 0, 0.0, 0.0 }
231 };
232 
233 #ifndef R2D
234 #  define R2D (180/M_PI)
235 #endif
236 #ifndef D2R
237 #  define D2R (M_PI/180)
238 #endif
239 
240 /* ==================================================================== */
241 /*      Some "standard" std::strings.                                        */
242 /* ==================================================================== */
243 
244 static const char ILW_False_Easting[] = "False Easting";
245 static const char ILW_False_Northing[] = "False Northing";
246 static const char ILW_Central_Meridian[] = "Central Meridian";
247 static const char ILW_Central_Parallel[] = "Central Parallel";
248 static const char ILW_Standard_Parallel_1[] = "Standard Parallel 1";
249 static const char ILW_Standard_Parallel_2[] = "Standard Parallel 2";
250 static const char ILW_Scale_Factor[] = "Scale Factor";
251 static const char ILW_Latitude_True_Scale[] = "Latitude of True Scale";
252 static const char ILW_Height_Persp_Center[] = "Height Persp. Center";
253 
ReadPrjParams(const std::string & section,const std::string & entry,const std::string & filename)254 static double ReadPrjParams(const std::string& section, const std::string& entry, const std::string& filename)
255 {
256     std::string str = ReadElement(section, entry, filename);
257     //string str="";
258     if (!str.empty())
259         return CPLAtof(str.c_str());
260 
261     return 0.0;
262 }
263 
fetchParams(const std::string & csyFileName,double * padfPrjParams)264 static int fetchParams(const std::string& csyFileName, double * padfPrjParams)
265 {
266     //Fill all projection parameters with zero
267     for ( int i = 0; i < 13; i++ )
268         padfPrjParams[i] = 0.0;
269 
270     //std::string pszProj = ReadElement("CoordSystem", "Projection", csyFileName);
271     std::string pszEllips = ReadElement("CoordSystem", "Ellipsoid", csyFileName);
272 
273     //fetch info about a custom ellipsoid
274     if( STARTS_WITH_CI(pszEllips.c_str(), "User Defined") )
275     {
276         padfPrjParams[0] = ReadPrjParams("Ellipsoid", "a", csyFileName);
277         padfPrjParams[2] = ReadPrjParams("Ellipsoid", "1/f", csyFileName);
278     }
279     else if( STARTS_WITH_CI(pszEllips.c_str(), "Sphere") )
280     {
281         padfPrjParams[0] = ReadPrjParams("CoordSystem", "Sphere Radius", csyFileName);
282     }
283 
284     padfPrjParams[3] = ReadPrjParams("Projection", "False Easting", csyFileName);
285     padfPrjParams[4] = ReadPrjParams("Projection", "False Northing", csyFileName);
286 
287     padfPrjParams[5] = ReadPrjParams("Projection", "Central Parallel", csyFileName);
288     padfPrjParams[6] = ReadPrjParams("Projection", "Central Meridian", csyFileName);
289 
290     padfPrjParams[7] = ReadPrjParams("Projection", "Standard Parallel 1", csyFileName);
291     padfPrjParams[8] = ReadPrjParams("Projection", "Standard Parallel 2", csyFileName);
292 
293     padfPrjParams[9] = ReadPrjParams("Projection", "Scale Factor", csyFileName);
294     padfPrjParams[10] = ReadPrjParams("Projection", "Latitude of True Scale", csyFileName);
295     padfPrjParams[11] = ReadPrjParams("Projection", "Zone", csyFileName);
296     padfPrjParams[12] = ReadPrjParams("Projection", ILW_Height_Persp_Center, csyFileName);
297 
298     return true;
299 }
300 
301 /************************************************************************/
302 /*                          mapTMParams                                  */
303 /************************************************************************/
304 /**
305  * fetch the parameters from ILWIS projection definition for
306  * --- Gauss-Krueger Germany.
307  * --- Gauss Colombia
308  * --- Gauss-Boaga Italy
309 **/
mapTMParams(std::string sProj,double dfZone,double & dfFalseEasting,double & dfCentralMeridian)310 static int mapTMParams(std::string sProj, double dfZone, double &dfFalseEasting, double &dfCentralMeridian)
311 {
312     if( STARTS_WITH_CI(sProj.c_str(), "Gauss-Krueger Germany") )
313     {
314         //Zone number must be in the range 1 to 3
315         dfCentralMeridian = 6.0 + (dfZone - 1) * 3;
316         dfFalseEasting = 2500000 + (dfZone - 1) * 1000000;
317     }
318     else if( STARTS_WITH_CI(sProj.c_str(), "Gauss-Boaga Italy") )
319     {
320         if ( dfZone == 1)
321         {
322             dfCentralMeridian = 9;
323             dfFalseEasting = 1500000;
324         }
325         else if ( dfZone == 2)
326         {
327             dfCentralMeridian = 15;
328             dfFalseEasting = 2520000;
329         }
330         else
331             return false;
332     }
333     else if( STARTS_WITH_CI(sProj.c_str(), "Gauss Colombia") )
334     {
335         //Zone number must be in the range 1 to 4
336         dfCentralMeridian = -77.08097220 + (dfZone - 1) * 3;
337     }
338     return true;
339 }
340 
341 /************************************************************************/
342 /*                          scaleFromLATTS()                             */
343 /************************************************************************/
344 /**
345  * Compute the scale factor from Latitude_Of_True_Scale parameter.
346  *
347 **/
scaleFromLATTS(std::string sEllips,double phits,double & scale)348 static void scaleFromLATTS( std::string sEllips, double phits, double &scale )
349 {
350     if( STARTS_WITH_CI(sEllips.c_str(), "Sphere") )
351     {
352         scale = cos(phits);
353     }
354     else
355     {
356         const IlwisEllips *piwEllips =  iwEllips;
357         double e2 = 0.0;
358         while ( piwEllips->pszIlwisEllips )
359         {
360             if( EQUALN( sEllips.c_str(), piwEllips->pszIlwisEllips, strlen(piwEllips->pszIlwisEllips) ) )
361             {
362                 double a = piwEllips->semiMajor;
363                 double b = a * ( 1 - piwEllips->invFlattening);
364                 e2 = ( a*a - b*b ) /( a*a );
365                 break;
366             }
367             piwEllips++;
368         }
369         scale = cos(phits) / sqrt (1. - e2 * sin(phits) * sin(phits));
370     }
371 }
372 
373 /************************************************************************/
374 /*                          ReadProjection()                           */
375 /************************************************************************/
376 
377 /**
378  * Import coordinate system from ILWIS projection definition.
379  *
380  * The method will import projection definition in ILWIS,
381  * It uses 13 parameters to define the coordinate system
382  * and datum/ellipsoid specified in the padfPrjParams array.
383  *
384  * @param csyFileName Name of .csy file
385 **/
386 
ReadProjection(const std::string & csyFileName)387 CPLErr ILWISDataset::ReadProjection( const std::string& csyFileName )
388 {
389     std::string pszEllips;
390     std::string pszDatum;
391     std::string pszProj;
392 
393     //translate ILWIS pre-defined coordinate systems
394     if( STARTS_WITH_CI(csyFileName.c_str(), "latlon.csy"))
395     {
396         pszProj = "LatLon";
397         pszDatum = "";
398         pszEllips = "Sphere";
399     }
400     else if ( STARTS_WITH_CI(csyFileName.c_str(), "LatlonWGS84.csy"))
401     {
402         pszProj = "LatLon";
403         pszDatum = "WGS 1984";
404         pszEllips = "WGS 84";
405     }
406     else
407     {
408         pszProj = ReadElement("CoordSystem", "Type", csyFileName);
409         if( !STARTS_WITH_CI( pszProj.c_str(), "LatLon" ) )
410             pszProj = ReadElement("CoordSystem", "Projection", csyFileName);
411         pszDatum = ReadElement("CoordSystem", "Datum", csyFileName);
412         pszEllips = ReadElement("CoordSystem", "Ellipsoid", csyFileName);
413     }
414 
415 /* -------------------------------------------------------------------- */
416 /*      Fetch array containing 13 coordinate system parameters          */
417 /* -------------------------------------------------------------------- */
418     double     padfPrjParams[13];
419     fetchParams(csyFileName, padfPrjParams);
420 
421     OGRSpatialReference oSRS;
422 /* -------------------------------------------------------------------- */
423 /*      Operate on the basis of the projection name.                    */
424 /* -------------------------------------------------------------------- */
425     if( STARTS_WITH_CI(pszProj.c_str(), "LatLon") )
426     {
427         //set datum later
428     }
429     else if( STARTS_WITH_CI(pszProj.c_str(), "Albers EqualArea Conic") )
430     {
431         oSRS.SetProjCS("Albers EqualArea Conic");
432         oSRS.SetACEA( padfPrjParams[7], padfPrjParams[8],
433                       padfPrjParams[5], padfPrjParams[6],
434                       padfPrjParams[3], padfPrjParams[4] );
435     }
436     else if( STARTS_WITH_CI(pszProj.c_str(), "Azimuthal Equidistant") )
437     {
438         oSRS.SetProjCS("Azimuthal Equidistant");
439         oSRS.SetAE( padfPrjParams[5], padfPrjParams[6],
440                     padfPrjParams[3], padfPrjParams[4] );
441     }
442     else if( STARTS_WITH_CI(pszProj.c_str(), "Central Cylindrical") )
443     {
444         //Use Central Parallel for dfStdP1
445         //padfPrjParams[5] is always to zero
446         oSRS.SetProjCS("Central Cylindrical");
447         oSRS.SetCEA( padfPrjParams[5], padfPrjParams[6],
448                      padfPrjParams[3], padfPrjParams[4] );
449     }
450     else if( STARTS_WITH_CI(pszProj.c_str(), "Cassini") )
451     {
452         //Use Latitude_Of_True_Scale for dfCenterLat
453         //Scale Factor 1.0 should always be defined
454         oSRS.SetProjCS("Cassini");
455         oSRS.SetCS(  padfPrjParams[10], padfPrjParams[6],
456                      padfPrjParams[3], padfPrjParams[4] );
457     }
458     else if( STARTS_WITH_CI(pszProj.c_str(), "DutchRD") )
459     {
460         oSRS.SetProjCS("DutchRD");
461         oSRS.SetStereographic  (  52.156160556,  5.387638889,
462                                   0.9999079,
463                                   155000,  463000);
464     }
465     else if( STARTS_WITH_CI(pszProj.c_str(), "Equidistant Conic") )
466     {
467         oSRS.SetProjCS("Equidistant Conic");
468         oSRS.SetEC(  padfPrjParams[7], padfPrjParams[8],
469                      padfPrjParams[5], padfPrjParams[6],
470                      padfPrjParams[3], padfPrjParams[4] );
471     }
472     else if( STARTS_WITH_CI(pszProj.c_str(), "Gauss-Krueger Germany") )
473     {
474         //FalseNorthing and CenterLat are always set to 0
475         //Scale 1.0 is defined
476         //FalseEasting and CentralMeridian are defined by the selected zone
477         mapTMParams("Gauss-Krueger Germany", padfPrjParams[11],
478                    padfPrjParams[3], padfPrjParams[6]);
479         oSRS.SetProjCS("Gauss-Krueger Germany");
480         oSRS.SetTM(  0, padfPrjParams[6],
481                      1.0,
482                      padfPrjParams[3], 0 );
483     }
484     else if ( STARTS_WITH_CI(pszProj.c_str(), "Gauss-Boaga Italy") )
485     {
486         //FalseNorthing and CenterLat are always set to 0
487         //Scale 0.9996 is defined
488         //FalseEasting and CentralMeridian are defined by the selected zone
489         mapTMParams("Gauss-Boaga Italy", padfPrjParams[11],
490                    padfPrjParams[3], padfPrjParams[6]);
491         oSRS.SetProjCS("Gauss-Boaga Italy");
492         oSRS.SetTM(  0, padfPrjParams[6],
493                      0.9996,
494                      padfPrjParams[3], 0 );
495     }
496     else if ( STARTS_WITH_CI(pszProj.c_str(), "Gauss Colombia"))
497     {
498         // 1000000 used for FalseNorthing and FalseEasting
499         // 1.0 used for scale
500         // CenterLat is defined 45.1609259259259
501         // CentralMeridian is defined by the selected zone
502         mapTMParams("Gauss Colombia", padfPrjParams[11],
503                    padfPrjParams[3], padfPrjParams[6]);
504         oSRS.SetProjCS("Gauss Colombia");
505         oSRS.SetTM(  45.1609259259259, padfPrjParams[6],
506                      1.0,
507                      1000000, 1000000 );
508     }
509     else if( STARTS_WITH_CI(pszProj.c_str(), "Gnomonic") )
510     {
511         oSRS.SetProjCS("Gnomonic");
512         oSRS.SetGnomonic( padfPrjParams[5], padfPrjParams[6],
513                           padfPrjParams[3], padfPrjParams[4] );
514     }
515     else if( STARTS_WITH_CI(pszProj.c_str(), "Lambert Conformal Conic") )
516     {
517         // should use 1.0 for scale factor in Ilwis definition
518         oSRS.SetProjCS("Lambert Conformal Conic");
519         oSRS.SetLCC( padfPrjParams[7], padfPrjParams[8],
520                      padfPrjParams[5], padfPrjParams[6],
521                      padfPrjParams[3], padfPrjParams[4] );
522     }
523     else if( STARTS_WITH_CI(pszProj.c_str(), "Lambert Cylind EqualArea") )
524     {
525         // Latitude_Of_True_Scale used for dfStdP1 ?
526         oSRS.SetProjCS("Lambert Conformal Conic");
527         oSRS.SetCEA( padfPrjParams[10],
528                      padfPrjParams[6],
529                      padfPrjParams[3], padfPrjParams[4] );
530     }
531     else if( STARTS_WITH_CI(pszProj.c_str(), "Mercator") )
532     {
533         // use 0 for CenterLat, scale is computed from the
534         // Latitude_Of_True_Scale
535         scaleFromLATTS( pszEllips, padfPrjParams[10], padfPrjParams[9] );
536         oSRS.SetProjCS("Mercator");
537         oSRS.SetMercator( 0, padfPrjParams[6],
538                           padfPrjParams[9],
539                           padfPrjParams[3], padfPrjParams[4] );
540     }
541     else if( STARTS_WITH_CI(pszProj.c_str(), "Miller") )
542     {
543         // use 0 for CenterLat
544         oSRS.SetProjCS("Miller");
545         oSRS.SetMC( 0, padfPrjParams[6],
546                     padfPrjParams[3], padfPrjParams[4] );
547     }
548     else if( STARTS_WITH_CI(pszProj.c_str(), "Mollweide") )
549     {
550         oSRS.SetProjCS("Mollweide");
551         oSRS.SetMollweide( padfPrjParams[6],
552                            padfPrjParams[3], padfPrjParams[4] );
553     }
554     else if( STARTS_WITH_CI(pszProj.c_str(), "Orthographic") )
555     {
556         oSRS.SetProjCS("Orthographic");
557         oSRS.SetOrthographic( padfPrjParams[5], padfPrjParams[6],
558                               padfPrjParams[3], padfPrjParams[4] );
559     }
560     else if( STARTS_WITH_CI(pszProj.c_str(), "Plate Carree") ||
561              STARTS_WITH_CI(pszProj.c_str(), "Plate Rectangle"))
562     {
563         // set 0.0 for CenterLat for Plate Carree projection
564         // skip Latitude_Of_True_Scale for Plate Rectangle projection definition
565         oSRS.SetProjCS(pszProj.c_str());
566         oSRS.SetEquirectangular( padfPrjParams[5], padfPrjParams[6],
567                                  padfPrjParams[3], padfPrjParams[4] );
568     }
569     else if( STARTS_WITH_CI(pszProj.c_str(), "PolyConic") )
570     {
571         // skip scale factor
572         oSRS.SetProjCS("PolyConic");
573         oSRS.SetPolyconic( padfPrjParams[5], padfPrjParams[6],
574                            padfPrjParams[3], padfPrjParams[4] );
575     }
576     else if( STARTS_WITH_CI(pszProj.c_str(), "Robinson") )
577     {
578         oSRS.SetProjCS("Robinson");
579         oSRS.SetRobinson( padfPrjParams[6],
580                           padfPrjParams[3], padfPrjParams[4] );
581     }
582     else if( STARTS_WITH_CI(pszProj.c_str(), "Sinusoidal") )
583     {
584         oSRS.SetProjCS("Sinusoidal");
585         oSRS.SetSinusoidal( padfPrjParams[6],
586                             padfPrjParams[3], padfPrjParams[4] );
587     }
588     else if( STARTS_WITH_CI(pszProj.c_str(), "Stereographic") )
589     {
590         oSRS.SetProjCS("Stereographic");
591         oSRS.SetStereographic( padfPrjParams[5], padfPrjParams[6],
592                                padfPrjParams[9],
593                                 padfPrjParams[3], padfPrjParams[4] );
594     }
595     else if( STARTS_WITH_CI(pszProj.c_str(), "Transverse Mercator") )
596     {
597         oSRS.SetProjCS("Transverse Mercator");
598         oSRS.SetStereographic( padfPrjParams[5], padfPrjParams[6],
599                                padfPrjParams[9],
600                                padfPrjParams[3], padfPrjParams[4] );
601     }
602     else if( STARTS_WITH_CI(pszProj.c_str(), "UTM") )
603     {
604         std::string pszNH = ReadElement("Projection", "Northern Hemisphere", csyFileName);
605         oSRS.SetProjCS("UTM");
606         if( STARTS_WITH_CI(pszNH.c_str(), "Yes") )
607             oSRS.SetUTM( (int) padfPrjParams[11], 1);
608         else
609             oSRS.SetUTM( (int) padfPrjParams[11], 0);
610     }
611     else if( STARTS_WITH_CI(pszProj.c_str(), "VanderGrinten") )
612     {
613         oSRS.SetVDG( padfPrjParams[6],
614                      padfPrjParams[3], padfPrjParams[4] );
615     }
616     else if( STARTS_WITH_CI(pszProj.c_str(), "GeoStationary Satellite") )
617     {
618         oSRS.SetGEOS( padfPrjParams[6],
619                       padfPrjParams[12],
620                       padfPrjParams[3],
621                       padfPrjParams[4] );
622     }
623     else if( STARTS_WITH_CI(pszProj.c_str(), "MSG Perspective") )
624     {
625         oSRS.SetGEOS( padfPrjParams[6],
626                       padfPrjParams[12],
627                       padfPrjParams[3],
628                       padfPrjParams[4] );
629     }
630     else
631     {
632         oSRS.SetLocalCS( pszProj.c_str() );
633     }
634 /* -------------------------------------------------------------------- */
635 /*      Try to translate the datum/spheroid.                            */
636 /* -------------------------------------------------------------------- */
637 
638     if ( !oSRS.IsLocal() )
639     {
640         const IlwisDatums   *piwDatum = iwDatums;
641 
642         // Search for matching datum
643         while ( piwDatum->pszIlwisDatum )
644         {
645             if( EQUALN( pszDatum.c_str(), piwDatum->pszIlwisDatum, strlen(piwDatum->pszIlwisDatum) ) )
646             {
647                 OGRSpatialReference oOGR;
648                 oOGR.importFromEPSG( piwDatum->nEPSGCode );
649                 oSRS.CopyGeogCSFrom( &oOGR );
650                 break;
651             }
652             piwDatum++;
653         } // End of searching for matching datum.
654 
655 /* -------------------------------------------------------------------- */
656 /*      If no matching for datum definition, fetch info about an        */
657 /*      ellipsoid.  semi major axis is always returned in meters        */
658 /* -------------------------------------------------------------------- */
659         const IlwisEllips *piwEllips =  iwEllips;
660         if (pszEllips.empty())
661             pszEllips="Sphere";
662         if ( !piwDatum->pszIlwisDatum )
663         {
664             while ( piwEllips->pszIlwisEllips )
665             {
666                 if( EQUALN( pszEllips.c_str(), piwEllips->pszIlwisEllips, strlen(piwEllips->pszIlwisEllips) ) )
667                 {
668                     double dfSemiMajor = piwEllips->semiMajor;
669                     if( STARTS_WITH_CI(pszEllips.c_str(), "Sphere") && padfPrjParams[0] != 0 )
670                     {
671                         dfSemiMajor = padfPrjParams[0];
672                     }
673                     oSRS.SetGeogCS( CPLSPrintf(
674                                         "Unknown datum based upon the %s ellipsoid",
675                                         piwEllips->pszIlwisEllips ),
676                                     CPLSPrintf(
677                                         "Not specified (based on %s spheroid)",
678                                         piwEllips->pszIlwisEllips ),
679                                     piwEllips->pszIlwisEllips,
680                                     dfSemiMajor,
681                                     piwEllips->invFlattening,
682                                     nullptr, 0.0, nullptr, 0.0 );
683                     oSRS.SetAuthority( "SPHEROID", "EPSG", piwEllips->nEPSGCode );
684 
685                     break;
686                 }
687                 piwEllips++;
688             } //end of searching for matching ellipsoid
689         }
690 
691 /* -------------------------------------------------------------------- */
692 /*      If no matching for ellipsoid definition, fetch info about an    */
693 /*      user defined ellipsoid. If cannot find, default to WGS 84       */
694 /* -------------------------------------------------------------------- */
695         if ( !piwEllips->pszIlwisEllips )
696         {
697 
698             if( STARTS_WITH_CI(pszEllips.c_str(), "User Defined") )
699             {
700                 oSRS.SetGeogCS( "Unknown datum based upon the custom ellipsoid",
701                                 "Not specified (based on custom ellipsoid)",
702                                 "Custom ellipsoid",
703                                 padfPrjParams[0], padfPrjParams[2],
704                                 nullptr, 0, nullptr, 0 );
705             }
706             else
707             {
708                 //if cannot find the user defined ellips, default to WGS84
709                 oSRS.SetWellKnownGeogCS( "WGS84" );
710             }
711         }
712     } // end of if ( !IsLocal() )
713 
714 /* -------------------------------------------------------------------- */
715 /*      Units translation                                          */
716 /* -------------------------------------------------------------------- */
717     if( oSRS.IsLocal() || oSRS.IsProjected() )
718     {
719         oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
720     }
721 
722     CPLFree(pszProjection);
723     oSRS.exportToWkt( &pszProjection );
724 
725     return CE_None;
726 }
727 
WriteFalseEastNorth(const std::string & csFileName,const OGRSpatialReference & oSRS)728 static void WriteFalseEastNorth(const std::string& csFileName, const OGRSpatialReference& oSRS)
729 {
730     WriteElement("Projection", ILW_False_Easting, csFileName,
731                  oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0));
732     WriteElement("Projection", ILW_False_Northing, csFileName,
733                  oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0));
734 }
735 
WriteProjectionName(const std::string & csFileName,const std::string & stProjection)736 static void WriteProjectionName(const std::string& csFileName, const std::string& stProjection)
737 {
738     WriteElement("CoordSystem", "Type", csFileName, "Projection");
739     WriteElement("CoordSystem", "Projection", csFileName, stProjection);
740 }
741 
WriteUTM(const std::string & csFileName,const OGRSpatialReference & oSRS)742 static void WriteUTM(const std::string& csFileName, const OGRSpatialReference& oSRS)
743 {
744     int bNorth;
745 
746     int nZone = oSRS.GetUTMZone( &bNorth );
747     WriteElement("CoordSystem", "Type", csFileName, "Projection");
748     WriteElement("CoordSystem", "Projection", csFileName, "UTM");
749     if (bNorth)
750         WriteElement("Projection", "Northern Hemisphere", csFileName, "Yes");
751     else
752         WriteElement("Projection", "Northern Hemisphere", csFileName, "No");
753     WriteElement("Projection", "Zone", csFileName, nZone);
754 }
755 
WriteAlbersConicEqualArea(const std::string & csFileName,const OGRSpatialReference & oSRS)756 static void WriteAlbersConicEqualArea(const std::string& csFileName, const OGRSpatialReference& oSRS)
757 {
758     WriteProjectionName(csFileName, "Albers EqualArea Conic");
759     WriteFalseEastNorth(csFileName, oSRS);
760     WriteElement("Projection", ILW_Central_Meridian, csFileName,
761                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
762     WriteElement("Projection", ILW_Central_Parallel, csFileName,
763                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
764     WriteElement("Projection", ILW_Standard_Parallel_1, csFileName,
765                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0));
766     WriteElement("Projection", ILW_Standard_Parallel_2, csFileName,
767                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0));
768 }
769 
WriteAzimuthalEquidistant(const std::string & csFileName,const OGRSpatialReference & oSRS)770 static void WriteAzimuthalEquidistant(const std::string& csFileName, const OGRSpatialReference& oSRS)
771 {
772     WriteProjectionName(csFileName, "Azimuthal Equidistant");
773     WriteFalseEastNorth(csFileName, oSRS);
774     WriteElement("Projection", ILW_Central_Meridian, csFileName,
775                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
776     WriteElement("Projection", ILW_Central_Parallel, csFileName,
777                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
778     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
779 }
780 
WriteCylindricalEqualArea(const std::string & csFileName,const OGRSpatialReference & oSRS)781 static void WriteCylindricalEqualArea(const std::string& csFileName, const OGRSpatialReference& oSRS)
782 {
783     WriteProjectionName(csFileName, "Central Cylindrical");
784     WriteFalseEastNorth(csFileName, oSRS);
785     WriteElement("Projection", ILW_Central_Meridian, csFileName,
786                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
787 }
788 
WriteCassiniSoldner(const std::string & csFileName,const OGRSpatialReference & oSRS)789 static void WriteCassiniSoldner(const std::string& csFileName, const OGRSpatialReference& oSRS)
790 {
791     WriteProjectionName(csFileName, "Cassini");
792     WriteFalseEastNorth(csFileName, oSRS);
793     WriteElement("Projection", ILW_Central_Meridian, csFileName,
794                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
795     WriteElement("Projection", ILW_Latitude_True_Scale, csFileName,
796                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
797     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
798 }
799 
WriteStereographic(const std::string & csFileName,const OGRSpatialReference & oSRS)800 static void WriteStereographic(const std::string& csFileName, const OGRSpatialReference& oSRS)
801 {
802     WriteProjectionName(csFileName, "Stereographic");
803     WriteFalseEastNorth(csFileName, oSRS);
804     WriteElement("Projection", ILW_Central_Meridian, csFileName,
805                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
806     WriteElement("Projection", ILW_Central_Parallel, csFileName,
807                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
808     WriteElement("Projection", ILW_Scale_Factor, csFileName,
809                  oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0));
810 }
811 
WriteEquidistantConic(const std::string & csFileName,const OGRSpatialReference & oSRS)812 static void WriteEquidistantConic(const std::string& csFileName, const OGRSpatialReference& oSRS)
813 {
814     WriteProjectionName(csFileName, "Equidistant Conic");
815     WriteFalseEastNorth(csFileName, oSRS);
816     WriteElement("Projection", ILW_Central_Meridian, csFileName,
817                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
818     WriteElement("Projection", ILW_Central_Parallel, csFileName,
819                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
820     WriteElement("Projection", ILW_Standard_Parallel_1, csFileName,
821                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0));
822     WriteElement("Projection", ILW_Standard_Parallel_2, csFileName,
823                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0));
824 }
825 
WriteTransverseMercator(const std::string & csFileName,const OGRSpatialReference & oSRS)826 static void WriteTransverseMercator(const std::string& csFileName, const OGRSpatialReference& oSRS)
827 {
828     WriteProjectionName(csFileName, "Transverse Mercator");
829     WriteFalseEastNorth(csFileName, oSRS);
830     WriteElement("Projection", ILW_Central_Meridian, csFileName,
831                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
832     WriteElement("Projection", ILW_Central_Parallel, csFileName,
833                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
834     WriteElement("Projection", ILW_Scale_Factor, csFileName,
835                  oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0));
836 }
837 
WriteGnomonic(const std::string & csFileName,const OGRSpatialReference & oSRS)838 static void WriteGnomonic(const std::string& csFileName, const OGRSpatialReference& oSRS)
839 {
840     WriteProjectionName(csFileName, "Gnomonic");
841     WriteFalseEastNorth(csFileName, oSRS);
842     WriteElement("Projection", ILW_Central_Meridian, csFileName,
843                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
844     WriteElement("Projection", ILW_Central_Parallel, csFileName,
845                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
846 }
847 
WriteLambertConformalConic(const std::string & csFileName,const OGRSpatialReference & oSRS)848 static void WriteLambertConformalConic(const std::string& csFileName, const OGRSpatialReference& oSRS)
849 {
850     WriteProjectionName(csFileName, "Lambert Conformal Conic");
851     WriteFalseEastNorth(csFileName, oSRS);
852     WriteElement("Projection", ILW_Central_Meridian, csFileName,
853                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
854     WriteElement("Projection", ILW_Central_Parallel, csFileName,
855                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
856     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
857 }
858 
WriteLambertConformalConic2SP(const std::string & csFileName,const OGRSpatialReference & oSRS)859 static void WriteLambertConformalConic2SP(const std::string& csFileName, const OGRSpatialReference& oSRS)
860 {
861     WriteProjectionName(csFileName, "Lambert Conformal Conic");
862     WriteFalseEastNorth(csFileName, oSRS);
863     WriteElement("Projection", ILW_Central_Meridian, csFileName,
864                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
865     WriteElement("Projection", ILW_Central_Parallel, csFileName,
866                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
867     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
868     WriteElement("Projection", ILW_Standard_Parallel_1, csFileName,
869                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0));
870     WriteElement("Projection", ILW_Standard_Parallel_2, csFileName,
871                  oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0));
872 }
873 
WriteLambertAzimuthalEqualArea(const std::string & csFileName,const OGRSpatialReference & oSRS)874 static void WriteLambertAzimuthalEqualArea(const std::string& csFileName, const OGRSpatialReference& oSRS)
875 {
876     WriteProjectionName(csFileName, "Lambert Azimuthal EqualArea");
877     WriteFalseEastNorth(csFileName, oSRS);
878     WriteElement("Projection", ILW_Central_Meridian, csFileName,
879                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
880     WriteElement("Projection", ILW_Central_Parallel, csFileName,
881                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
882 }
883 
WriteMercator_1SP(const std::string & csFileName,const OGRSpatialReference & oSRS)884 static void WriteMercator_1SP(const std::string& csFileName, const OGRSpatialReference& oSRS)
885 {
886     WriteProjectionName(csFileName, "Mercator");
887     WriteFalseEastNorth(csFileName, oSRS);
888     WriteElement("Projection", ILW_Central_Meridian, csFileName,
889                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
890     WriteElement("Projection", ILW_Latitude_True_Scale, csFileName,
891                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
892 }
893 
WriteMillerCylindrical(const std::string & csFileName,const OGRSpatialReference & oSRS)894 static void WriteMillerCylindrical(const std::string& csFileName, const OGRSpatialReference& oSRS)
895 {
896     WriteProjectionName(csFileName, "Miller");
897     WriteFalseEastNorth(csFileName, oSRS);
898     WriteElement("Projection", ILW_Central_Meridian, csFileName,
899                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
900 }
901 
WriteMolleweide(const std::string & csFileName,const OGRSpatialReference & oSRS)902 static void WriteMolleweide(const std::string& csFileName, const OGRSpatialReference& oSRS)
903 {
904     WriteProjectionName(csFileName, "Mollweide");
905     WriteFalseEastNorth(csFileName, oSRS);
906     WriteElement("Projection", ILW_Central_Meridian, csFileName,
907                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
908 }
909 
WriteOrthographic(const std::string & csFileName,const OGRSpatialReference & oSRS)910 static void WriteOrthographic(const std::string& csFileName, const OGRSpatialReference& oSRS)
911 {
912     WriteProjectionName(csFileName, "Orthographic");
913     WriteFalseEastNorth(csFileName, oSRS);
914     WriteElement("Projection", ILW_Central_Meridian, csFileName,
915                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
916     WriteElement("Projection", ILW_Central_Parallel, csFileName,
917                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
918 }
919 
WritePlateRectangle(const std::string & csFileName,const OGRSpatialReference & oSRS)920 static void WritePlateRectangle(const std::string& csFileName, const OGRSpatialReference& oSRS)
921 {
922     WriteProjectionName(csFileName, "Plate Rectangle");
923     WriteFalseEastNorth(csFileName, oSRS);
924     WriteElement("Projection", ILW_Central_Meridian, csFileName,
925                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
926     WriteElement("Projection", ILW_Central_Parallel, csFileName,
927                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
928     WriteElement("Projection", ILW_Latitude_True_Scale, csFileName, "0.0000000000");
929 }
930 
WritePolyConic(const std::string & csFileName,const OGRSpatialReference & oSRS)931 static void WritePolyConic(const std::string& csFileName, const OGRSpatialReference& oSRS)
932 {
933     WriteProjectionName(csFileName, "PolyConic");
934     WriteFalseEastNorth(csFileName, oSRS);
935     WriteElement("Projection", ILW_Central_Meridian, csFileName,
936                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
937     WriteElement("Projection", ILW_Central_Parallel, csFileName,
938                  oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
939     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
940 }
941 
WriteRobinson(const std::string & csFileName,const OGRSpatialReference & oSRS)942 static void WriteRobinson(const std::string& csFileName, const OGRSpatialReference& oSRS)
943 {
944     WriteProjectionName(csFileName, "Robinson");
945     WriteFalseEastNorth(csFileName, oSRS);
946     WriteElement("Projection", ILW_Central_Meridian, csFileName,
947                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
948 }
949 
WriteSinusoidal(const std::string & csFileName,const OGRSpatialReference & oSRS)950 static void WriteSinusoidal(const std::string& csFileName, const OGRSpatialReference& oSRS)
951 {
952     WriteProjectionName(csFileName, "Sinusoidal");
953     WriteFalseEastNorth(csFileName, oSRS);
954     WriteElement("Projection", ILW_Central_Meridian, csFileName,
955                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
956 }
957 
WriteVanderGrinten(const std::string & csFileName,const OGRSpatialReference & oSRS)958 static void WriteVanderGrinten(const std::string& csFileName, const OGRSpatialReference& oSRS)
959 {
960     WriteProjectionName(csFileName, "VanderGrinten");
961     WriteFalseEastNorth(csFileName, oSRS);
962     WriteElement("Projection", ILW_Central_Meridian, csFileName,
963                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
964 }
965 
WriteGeoStatSat(const std::string & csFileName,const OGRSpatialReference & oSRS)966 static void WriteGeoStatSat(const std::string& csFileName, const OGRSpatialReference& oSRS)
967 {
968     WriteProjectionName(csFileName, "GeoStationary Satellite");
969     WriteFalseEastNorth(csFileName, oSRS);
970     WriteElement("Projection", ILW_Central_Meridian, csFileName,
971                  oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0));
972     WriteElement("Projection", ILW_Scale_Factor, csFileName, "1.0000000000");
973     WriteElement("Projection", ILW_Height_Persp_Center, csFileName,
974                  oSRS.GetNormProjParm(SRS_PP_SATELLITE_HEIGHT, 35785831.0));
975 }
976 
977 /************************************************************************/
978 /*                          WriteProjection()                           */
979 /************************************************************************/
980 /**
981  * Export coordinate system in ILWIS projection definition.
982  *
983  * Converts the loaded coordinate reference system into ILWIS projection
984  * definition to the extent possible.
985  */
WriteProjection()986 CPLErr ILWISDataset::WriteProjection()
987 
988 {
989     OGRSpatialReference oSRS;
990     OGRSpatialReference *poGeogSRS = nullptr;
991 
992     std::string csFileName = CPLResetExtension(osFileName, "csy" );
993     std::string pszBaseName = std::string(CPLGetBasename( osFileName ));
994     //std::string pszPath = std::string(CPLGetPath( osFileName ));
995     bool bProjection = ((pszProjection != nullptr) && (strlen(pszProjection)>0));
996     bool bHaveSRS;
997     if( bProjection && (oSRS.importFromWkt( pszProjection ) == OGRERR_NONE) )
998     {
999         bHaveSRS = true;
1000     }
1001     else
1002         bHaveSRS = false;
1003 
1004     const IlwisDatums   *piwDatum = iwDatums;
1005     //std::string pszEllips;
1006     std::string osDatum;
1007     //std::string pszProj;
1008 
1009 /* -------------------------------------------------------------------- */
1010 /*      Collect datum/ellips information.                                      */
1011 /* -------------------------------------------------------------------- */
1012     if( bHaveSRS )
1013     {
1014         poGeogSRS = oSRS.CloneGeogCS();
1015     }
1016 
1017     std::string grFileName = CPLResetExtension(osFileName, "grf" );
1018     std::string csy;
1019     if( poGeogSRS )
1020     {
1021         csy = pszBaseName + ".csy";
1022 
1023         WriteElement("Ilwis", "Type", csFileName, "CoordSystem");
1024         const char* pszDatum = poGeogSRS->GetAttrValue( "GEOGCS|DATUM" );
1025         if( pszDatum )
1026             osDatum = pszDatum;
1027 
1028         /* WKT to ILWIS translation */
1029         while ( piwDatum->pszWKTDatum)
1030         {
1031             if( EQUALN( osDatum.c_str(), piwDatum->pszWKTDatum, strlen(piwDatum->pszWKTDatum) ) )
1032             {
1033                 WriteElement("CoordSystem", "Datum", csFileName, piwDatum->pszIlwisDatum);
1034                 break;
1035             }
1036             piwDatum++;
1037         } // End of searching for matching datum.
1038         WriteElement("CoordSystem", "Width", csFileName, 28);
1039        // pszEllips = poGeogSRS->GetAttrValue( "GEOGCS|DATUM|SPHEROID" );
1040         double a = poGeogSRS->GetSemiMajor();
1041         /* b = */ poGeogSRS->GetSemiMinor();
1042         double f = poGeogSRS->GetInvFlattening();
1043         WriteElement("CoordSystem", "Ellipsoid", csFileName, "User Defined");
1044         WriteElement("Ellipsoid", "a", csFileName, a);
1045         WriteElement("Ellipsoid", "1/f", csFileName, f);
1046     }
1047     else
1048         csy = "unknown.csy";
1049 
1050 /* -------------------------------------------------------------------- */
1051 /*  Determine to write a geo-referencing file for the dataset to create */
1052 /* -------------------------------------------------------------------- */
1053     if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
1054         || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
1055         || adfGeoTransform[4] != 0.0 || fabs(adfGeoTransform[5]) != 1.0)
1056         WriteElement("GeoRef", "CoordSystem", grFileName, csy);
1057 
1058 /* -------------------------------------------------------------------- */
1059 /*  Recognise various projections.                                      */
1060 /* -------------------------------------------------------------------- */
1061     const char * pszProjName = nullptr;
1062 
1063     if( bHaveSRS )
1064         pszProjName = oSRS.GetAttrValue( "PROJCS|PROJECTION" );
1065 
1066     if( pszProjName == nullptr )
1067     {
1068         if( bHaveSRS && oSRS.IsGeographic() )
1069         {
1070             WriteElement("CoordSystem", "Type", csFileName, "LatLon");
1071         }
1072     }
1073     else if( oSRS.GetUTMZone( nullptr ) != 0 )
1074     {
1075         WriteUTM(csFileName, oSRS);
1076     }
1077     else if( EQUAL(pszProjName,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
1078     {
1079         WriteAlbersConicEqualArea(csFileName, oSRS);
1080     }
1081     else if( EQUAL(pszProjName,SRS_PT_AZIMUTHAL_EQUIDISTANT) )
1082     {
1083         WriteAzimuthalEquidistant(csFileName, oSRS);
1084     }
1085     else if( EQUAL(pszProjName,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
1086     {
1087         WriteCylindricalEqualArea(csFileName, oSRS);
1088     }
1089     else if( EQUAL(pszProjName,SRS_PT_CASSINI_SOLDNER) )
1090     {
1091         WriteCassiniSoldner(csFileName, oSRS);
1092     }
1093     else if( EQUAL(pszProjName,SRS_PT_STEREOGRAPHIC) )
1094     {
1095         WriteStereographic(csFileName, oSRS);
1096     }
1097     else if( EQUAL(pszProjName,SRS_PT_EQUIDISTANT_CONIC) )
1098     {
1099         WriteEquidistantConic(csFileName, oSRS);
1100     }
1101     else if( EQUAL(pszProjName,SRS_PT_TRANSVERSE_MERCATOR) )
1102     {
1103         WriteTransverseMercator(csFileName, oSRS);
1104     }
1105     else if( EQUAL(pszProjName,SRS_PT_GNOMONIC) )
1106     {
1107         WriteGnomonic(csFileName, oSRS);
1108     }
1109     else if( EQUAL(pszProjName,"Lambert_Conformal_Conic") )
1110     {
1111         WriteLambertConformalConic(csFileName, oSRS);
1112     }
1113     else if( EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) )
1114     {
1115         WriteLambertConformalConic(csFileName, oSRS);
1116     }
1117     else if( EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
1118     {
1119         WriteLambertConformalConic2SP(csFileName, oSRS);
1120     }
1121     else if( EQUAL(pszProjName,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
1122     {
1123         WriteLambertAzimuthalEqualArea(csFileName, oSRS);
1124     }
1125     else if( EQUAL(pszProjName,SRS_PT_MERCATOR_1SP) )
1126     {
1127         WriteMercator_1SP(csFileName, oSRS);
1128     }
1129     else if( EQUAL(pszProjName,SRS_PT_MILLER_CYLINDRICAL) )
1130     {
1131         WriteMillerCylindrical(csFileName, oSRS);
1132     }
1133     else if( EQUAL(pszProjName,SRS_PT_MOLLWEIDE) )
1134     {
1135         WriteMolleweide(csFileName, oSRS);
1136     }
1137     else if( EQUAL(pszProjName,SRS_PT_ORTHOGRAPHIC) )
1138     {
1139         WriteOrthographic(csFileName, oSRS);
1140     }
1141     else if( EQUAL(pszProjName,SRS_PT_EQUIRECTANGULAR) )
1142     {
1143         WritePlateRectangle(csFileName, oSRS);
1144     }
1145     else if( EQUAL(pszProjName,SRS_PT_POLYCONIC) )
1146     {
1147         WritePolyConic(csFileName, oSRS);
1148     }
1149     else if( EQUAL(pszProjName,SRS_PT_ROBINSON) )
1150     {
1151         WriteRobinson(csFileName, oSRS);
1152     }
1153     else if( EQUAL(pszProjName,SRS_PT_SINUSOIDAL) )
1154     {
1155         WriteSinusoidal(csFileName, oSRS);
1156     }
1157     else if( EQUAL(pszProjName,SRS_PT_VANDERGRINTEN) )
1158     {
1159         WriteVanderGrinten(csFileName, oSRS);
1160     }
1161     else if( EQUAL(pszProjName,SRS_PT_GEOSTATIONARY_SATELLITE) )
1162     {
1163         WriteGeoStatSat(csFileName, oSRS);
1164     }
1165     else
1166     {
1167         // Projection unknown by ILWIS
1168     }
1169 
1170 /* -------------------------------------------------------------------- */
1171 /*      Cleanup                                                         */
1172 /* -------------------------------------------------------------------- */
1173     if( poGeogSRS != nullptr )
1174         delete poGeogSRS;
1175 
1176     return CE_None;
1177 }
1178 
1179 } // namespace GDAL
1180