Implement RFC 2: Initial integration of "GDAL SRS barn" work

This work mostly consists of:
- a C++ implementation of the ISO-19111:2018 / OGC Topic 2
  "Referencing by coordinates" classes to represent Datums,
  Coordinate systems, CRSs (Coordinate Reference Systems) and
  Coordinate Operations.
- methods to convert between this C++ modeling and WKT1, WKT2
  and PROJ string representations of those objects
- management and query of a SQLite3 database of CRS and Coordinate Operation definition
- a C API binding part of those capabilities

This is all-in-one squashed commit of the work of
https://github.com/OSGeo/proj.4/pull/1040

1 files changed, 527 insertions, 0 deletions

diff --git a/scripts/build_db_create_ignf.py b/scripts/build_db_create_ignf.py
new file mode 100755
index 00000000..97157ec5
--- /dev/null
+++ b/scripts/build_db_create_ignf.py
@@ -0,0 +1,527 @@
+#!/usr/bin/env python
+###############################################################################
+# $Id$
+#
+#  Project:  PROJ
+#  Purpose:  Build SRS and coordinate transform database
+#  Author:   Even Rouault <even.rouault at spatialys.com>
+#
+###############################################################################
+#  Copyright (c) 2018, Even Rouault <even.rouault at spatialys.com>
+#
+#  Permission is hereby granted, free of charge, to any person obtaining a
+#  copy of this software and associated documentation files (the "Software"),
+#  to deal in the Software without restriction, including without limitation
+#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+#  and/or sell copies of the Software, and to permit persons to whom the
+#  Software is furnished to do so, subject to the following conditions:
+#
+#  The above copyright notice and this permission notice shall be included
+#  in all copies or substantial portions of the Software.
+#
+#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+#  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+#  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+#  DEALINGS IN THE SOFTWARE.
+###############################################################################
+
+import copy
+import json
+import os
+import sqlite3
+import sys
+
+if len(sys.argv) != 3:
+    print('Usage: build_db_create_ignf.py path_to_IGNF_file proj.db')
+    sys.exit(1)
+
+def escape_literal(x):
+    return x.replace("'", "''")
+
+def find_ellipsoid(cursor, line, d):
+    for prec in (10,9,8,7,6):
+        cursor.execute('SELECT code, name FROM ellipsoid WHERE code != 1024 AND semi_major_axis = ? AND (inv_flattening = ? OR (abs(semi_major_axis / (semi_major_axis - semi_minor_axis)) - ?) < 1e-%d)' % prec, (float(d['a']), float(d['rf']), float(d['rf'])))
+        row = cursor.fetchone()
+        if prec > 6 and row is None:
+            continue
+        assert row, line
+        ellps_code, ellps_name = row
+        return ellps_code, ellps_name
+
+def find_geogcrs_code(line,d):
+    pm = None
+    if 'pm' in d:
+        pm = d['pm']
+        assert pm == '2.337229167', pm
+
+    key = { 'a' : d['a'], 'rf' : d['rf'], 'towgs84' : d['towgs84'], 'pm': pm }
+    key_geog_crs = json.dumps(key)
+    if 'W84' in code:
+        geogcrs_code = "WGS84G"
+    elif 'RGFG95' in code:
+        geogcrs_code = "RGFG95GEO"
+    elif 'ETRS89' in code:
+        geogcrs_code = "ETRS89GEO"
+    elif 'RGF93' in code:
+        geogcrs_code = "RGF93G"
+    elif 'LAMB93' in code:
+        geogcrs_code = "RGF93G"
+    elif 'MILLER' in code:
+        geogcrs_code = "WGS84G"
+    elif code.startswith('GEOPORTAL'):
+        geogcrs_code = 'RGF93G'
+    elif key_geog_crs not in map_geog_crs:
+        max_count_common = 0
+        max_common_code = None
+        for k in map_geog_crs:
+            temp = copy.copy(key)
+            temp['code'] = map_geog_crs[k]['code']
+            if temp == json.loads(k):
+                count_common = 0
+                while True:
+                    if len(code) <= count_common:
+                        break
+                    if len(map_geog_crs[k]['code']) <= count_common:
+                        break
+                    if code[count_common] == map_geog_crs[k]['code'][count_common]:
+                        count_common += 1
+                    else:
+                        break
+                if count_common > max_count_common:
+                    max_count_common = count_common
+                    max_common_code = map_geog_crs[k]['code']
+        assert max_count_common >= 4, (max_common_code, max_count_common, line)
+        geogcrs_code = max_common_code
+    else:
+        geogcrs_code = map_geog_crs[key_geog_crs]['code']
+    return geogcrs_code
+
+
+IGNF_file = sys.argv[1]
+proj_db = sys.argv[2]
+
+conn = sqlite3.connect(proj_db)
+cursor = conn.cursor()
+
+all_sql = []
+
+geocentric_crs_name_to_datum_code = {}
+map_geog_crs = {}
+
+for line in open(IGNF_file, 'rt').readlines():
+    line = line.strip('\n')
+    if line[0] == '#' or line.startswith('<metadata>'):
+        continue
+    assert line[0] == '<'
+    code = line[1:line.find('>')]
+    proj4_string = line[line.find('>')+1:]
+    #print(code, proj4_string)
+    tokens = proj4_string.split(' +')
+    assert tokens[0] == ''
+    tokens = tokens[1:]
+    d = {}
+    for token in tokens:
+        if token == 'no_defs <>' or token == 'wktext':
+            continue
+        pos_equal = token.find('=')
+        assert pos_equal > 0, line
+        key = token[0:pos_equal]
+        value = token[pos_equal+1:]
+        d[key] = value
+    assert 'proj' in d, line
+    assert 'title' in d, line
+    proj = d['proj']
+
+    assert 'a' in d
+    assert 'rf' in d
+    assert (not 'nadgrids' in d) or d['nadgrids'] in ('null', '@null', 'ntf_r93.gsb,null'), line
+    assert 'towgs84' in d
+
+    if proj == 'geocent':
+
+        ellps_code, ellps_name = find_ellipsoid(cursor, line, d)
+
+        datum_code = 'DATUM_' + code
+
+        geocentric_crs_name_to_datum_code[d['title']] = datum_code
+
+        sql = """INSERT INTO "geodetic_datum" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','%s','EPSG','8901','EPSG','1262',0);""" % (datum_code, escape_literal(d['title']), ellps_code)
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','geocentric');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "geodetic_crs" VALUES('IGNF','%s','%s',NULL,NULL,'geocentric','EPSG','6500','IGNF','%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), datum_code)
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4978','helmert_transformation');""" % (code)
+        #all_sql.append(sql)
+
+        towgs84 = d['towgs84'].split(',')
+        assert len(towgs84) in (3, 7), towgs84
+        x = towgs84[0]
+        y = towgs84[1]
+        z = towgs84[2]
+        if len(towgs84) == 3:
+            rx = 'NULL'
+            ry = 'NULL'
+            rz = 'NULL'
+            s = 'NULL'
+            r_uom_auth_name = 'NULL'
+            r_uom_code = 'NULL'
+            s_uom_auth_name = 'NULL'
+            s_uom_code = 'NULL'
+            method_code = "'1031'"
+            method_name = "'Geocentric translations (geocentric domain)'"
+        else:
+            rx = towgs84[3]
+            ry = towgs84[4]
+            rz = towgs84[5]
+            s = towgs84[6]
+            r_uom_auth_name = "'EPSG'"
+            r_uom_code = "'9104'"
+            s_uom_auth_name = "'EPSG'"
+            s_uom_code = "'9202'"
+            method_code = "'1033'"
+            method_name = "'Position Vector transformation (geocentric domain)'"
+
+        sql = """INSERT INTO "helmert_transformation" VALUES('IGNF','IGNF_%s_TO_EPSG_4978','%s to WGS 84',NULL,NULL,'EPSG',%s,%s,'IGNF','%s','EPSG','4978','EPSG','1262',NULL,%s,%s,%s,'EPSG','9001',%s,%s,%s,%s,%s,%s,%s, %s,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), method_code, method_name, code, x, y, z, rx, ry, rz, r_uom_auth_name, r_uom_code, s, s_uom_auth_name, s_uom_code)
+        all_sql.append(sql)
+
+        if 'nadgrids' in d and d['nadgrids'] == 'ntf_r93.gsb,null':
+
+            #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4978_GRID','grid_transformation');""" % (code)
+            #all_sql.append(sql)
+
+            sql = """INSERT INTO "grid_transformation" VALUES('IGNF','IGNF_%s_TO_EPSG_4978_GRID','%s to WGS 84 (2)',NULL,NULL,'EPSG','9615','NTv2','IGNF','%s','EPSG','4978','EPSG','1262',NULL,'EPSG','8656','Latitude and longitude difference file','ntf_r93.gsb',NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), code)
+            all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'longlat':
+
+        pm = None
+        pm_code = "'8901'" # Greenwich
+        if 'pm' in d:
+            pm = d['pm']
+            assert pm == '2.337229167', pm
+            pm_code = "'8903'" # Paris
+
+        ellps_code, ellps_name = find_ellipsoid(cursor, line, d)
+
+        if d['title'] in geocentric_crs_name_to_datum_code:
+            datum_code = geocentric_crs_name_to_datum_code[d['title']]
+        else:
+            datum_code = 'DATUM_' + code
+            sql = """INSERT INTO "geodetic_datum" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','%s','EPSG',%s,'EPSG','1262',0);""" % (datum_code, escape_literal(d['title']), ellps_code, pm_code)
+            all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','geographic 2D');""" % (code)
+        #all_sql.append(sql)
+
+        cs_code = "'6424'" # Long Lat deg
+        if 'grades' in d['title']:
+            cs_code = "'6425'" # Long Lat grad
+
+        sql = """INSERT INTO "geodetic_crs" VALUES('IGNF','%s','%s',NULL,NULL,'geographic 2D','EPSG',%s,'IGNF','%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), cs_code, datum_code)
+        all_sql.append(sql)
+
+        key = { 'a' : d['a'], 'rf' : d['rf'], 'towgs84' : d['towgs84'], 'pm': pm }
+        if "+towgs84=0.0000,0.0000,0.0000 +a=6378137.0000 +rf=298.2572221010000" in line or code == 'WGS84G':
+            key['code'] = code
+        key_geog_crs = json.dumps(key)
+        assert key_geog_crs not in map_geog_crs, (line, map_geog_crs[key_geog_crs])
+        map_geog_crs[ key_geog_crs ] = { 'code': code, 'title': d['title'] }
+
+        if code == 'NTFP': # Grades
+            assert cs_code == "'6425'" # Long Lat grad
+
+            #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_NTFP_TO_IGNF_NTFG','other_transformation');"""
+            #all_sql.append(sql)
+
+            sql = """INSERT INTO "other_transformation" VALUES('IGNF','IGNF_NTFP_TO_IGNF_NTFG','Nouvelle Triangulation Francaise Paris grades to Nouvelle Triangulation Francaise Greenwich degres sexagesimaux',NULL,NULL,'EPSG','9601','Longitude rotation','IGNF','NTFP','IGNF','NTFG','EPSG','1262',0.0,'EPSG','8602','Longitude offset',2.5969213,'EPSG','9105',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);"""
+            all_sql.append(sql)
+
+            #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326','concatenated_operation');""" % (code)
+            #all_sql.append(sql)
+
+            sql = """INSERT INTO "concatenated_operation" VALUES('IGNF','IGNF_NTFP_TO_EPSG_4326','Nouvelle Triangulation Francaise Paris grades to WGS 84',NULL,NULL,'IGNF','NTFP','EPSG','4326','EPSG','1262',NULL,'IGNF','IGNF_NTFP_TO_IGNF_NTFG','IGNF','IGNF_NTFG_TO_EPSG_4326',NULL,NULL,0);"""
+            all_sql.append(sql)
+
+            #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326_GRID','concatenated_operation');""" % (code)
+            #all_sql.append(sql)
+
+            sql = """INSERT INTO "concatenated_operation" VALUES('IGNF','IGNF_NTFP_TO_EPSG_4326_GRID','Nouvelle Triangulation Francaise Paris grades to WGS 84 (2)',NULL,NULL,'IGNF','NTFP','EPSG','4326','EPSG','1262',NULL,'IGNF','IGNF_NTFP_TO_IGNF_NTFG','IGNF','IGNF_NTFG_TO_EPSG_4326_GRID',NULL,NULL,0);"""
+            all_sql.append(sql)
+
+            continue
+
+        assert cs_code ==  "'6424'" # Long Lat deg
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326','helmert_transformation');""" % (code)
+        #all_sql.append(sql)
+
+        towgs84 = d['towgs84'].split(',')
+        assert len(towgs84) in (3, 7), towgs84
+        x = towgs84[0]
+        y = towgs84[1]
+        z = towgs84[2]
+        if len(towgs84) == 3:
+            rx = 'NULL'
+            ry = 'NULL'
+            rz = 'NULL'
+            s = 'NULL'
+            r_uom_auth_name = 'NULL'
+            r_uom_code = 'NULL'
+            s_uom_auth_name = 'NULL'
+            s_uom_code = 'NULL'
+            method_code = "'9603'"
+            method_name = "'Geocentric translations (geog2D domain)'"
+ 
+        else:
+            rx = towgs84[3]
+            ry = towgs84[4]
+            rz = towgs84[5]
+            s = towgs84[6]
+            r_uom_auth_name = "'EPSG'"
+            r_uom_code = "'9104'"
+            s_uom_auth_name = "'EPSG'"
+            s_uom_code = "'9202'"
+            method_code = "'9606'"
+            method_name = "'Position Vector transformation (geog2D domain)'"
+ 
+        sql = """INSERT INTO "helmert_transformation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326','%s to WGS 84',NULL,NULL,'EPSG',%s,%s,'IGNF','%s','EPSG','4326','EPSG','1262',NULL,%s,%s,%s,'EPSG','9001',%s,%s,%s,%s,%s,%s,%s, %s,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), method_code, method_name, code, x, y, z, rx, ry, rz, r_uom_auth_name, r_uom_code, s, s_uom_auth_name, s_uom_code)
+        all_sql.append(sql)
+
+        if 'nadgrids' in d and d['nadgrids'] == 'ntf_r93.gsb,null':
+
+            #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326_GRID','grid_transformation');""" % (code)
+            #all_sql.append(sql)
+
+            sql = """INSERT INTO "grid_transformation" VALUES('IGNF','IGNF_%s_TO_EPSG_4326_GRID','%s to WGS 84 (2)',NULL,NULL,'EPSG','9615','NTv2','IGNF','%s','EPSG','4326','EPSG','1262',NULL,'EPSG','8656','Latitude and longitude difference file','ntf_r93.gsb',NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), code)
+            all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'tmerc':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        assert 'k_0' in d
+        assert 'lat_0' in d
+        assert 'lon_0' in d
+        assert 'units' in d and d['units'] == 'm'
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',%s,'EPSG','9102','EPSG','8802','Longitude of natural origin',%s,'EPSG','9102','EPSG','8805','Scale factor at natural origin',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_0'], d['lon_0'], d['k_0'], d['x_0'], d['y_0'])
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'laea':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        assert 'lat_0' in d
+        assert 'lon_0' in d
+        assert 'units' in d and d['units'] == 'm'
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','9820','Lambert Azimuthal Equal Area','EPSG','8801','Latitude of natural origin',%s,'EPSG','9102','EPSG','8802','Longitude of natural origin',%s,'EPSG','9102','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_0'], d['lon_0'], d['x_0'], d['y_0'])
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'lcc':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        assert 'lat_0' in d
+        assert 'lat_1' in d
+        assert 'k_0' in d or 'lat_2' in d, line
+        assert 'lon_0' in d
+        assert 'units' in d and d['units'] == 'm'
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        if 'lat_2' in d:
+            sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','9802','Lambert Conic Conformal (2SP)','EPSG','8821','Latitude of false origin',%s,'EPSG','9102','EPSG','8822','Longitude of false origin',%s,'EPSG','9102','EPSG','8823','Latitude of 1st standard parallel',%s,'EPSG','9102','EPSG','8824','Latitude of 2nd standard parallel',%s,'EPSG','9102','EPSG','8826','Easting at false origin',%s,'EPSG','9001','EPSG','8827','Northing at false origin',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_0'], d['lon_0'], d['lat_1'], d['lat_2'], d['x_0'], d['y_0'])
+        else:
+            assert d['lat_0'] == d['lat_1'], line
+            sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','9801','Lambert Conic Conformal (1SP)','EPSG','8801','Latitude of natural origin',%s,'EPSG','9102','EPSG','8802','Longitude of natural origin',%s,'EPSG','9102','EPSG','8805','Scale factor at natural origin',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_0'], d['lon_0'], d['k_0'], d['x_0'], d['y_0'])
+
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'eqc':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d, line
+        assert 'y_0' in d, line
+        assert 'lat_0' in d, line
+        assert 'lat_ts' in d, line
+        assert 'lon_0' in d, line
+        assert 'units' in d and d['units'] == 'm', line
+        assert float(d['lat_0']) == 0, line
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','1028','Equidistant Cylindrical','EPSG','8823','Latitude of 1st standard parallel',%s,'EPSG','9102','EPSG','8802','Longitude of natural origin',%s,'EPSG','9102','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_ts'], d['lon_0'], d['x_0'], d['y_0'])
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'gstmerc':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        assert 'k_0' in d
+        assert 'lat_0' in d
+        assert 'lon_0' in d
+        assert 'units' in d and d['units'] == 'm'
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262',NULL,NULL,'Gauss Schreiber Transverse Mercator','EPSG','8801','Latitude of natural origin',%s,'EPSG','9102','EPSG','8802','Longitude of natural origin',%s,'EPSG','9102','EPSG','8805','Scale factor at natural origin',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_0'], d['lon_0'], d['k_0'], d['x_0'], d['y_0'])
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'stere':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        if 'k' in d or 'k_0' in d:
+            print('k/k_0 ignored for ' + line)
+        assert 'lat_0' in d
+        assert float(d['lat_0']) == -90, line
+        assert 'lat_ts' in d
+        assert 'lon_0' in d
+        assert float(d['lon_0']) == 140, line
+        assert 'units' in d and d['units'] == 'm'
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262','EPSG','9829','Polar Stereographic (variant B)','EPSG','8832','Latitude of standard parallel',%s,'EPSG','9102','EPSG','8833','Longitude of origin',%s,'EPSG','9102','EPSG','8806','False easting',%s,'EPSG','9001','EPSG','8807','False northing',%s,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']), d['lat_ts'], d['lon_0'], d['x_0'], d['y_0'])
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        # Hardcoding 1025 for coordinate system since the only instance uses it
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','1025','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    elif proj == 'mill':
+
+        geogcrs_code = find_geogcrs_code(line,d)
+
+        assert 'x_0' in d
+        assert 'y_0' in d
+        assert 'lon_0' in d
+        assert 'units' in d and d['units'] == 'm'
+        assert float(d['x_0']) == 0, d
+        assert float(d['y_0']) == 0, d
+        assert float(d['lon_0']) == 0, d
+
+        #sql = """INSERT INTO "coordinate_operation" VALUES('IGNF','CONV_%s','conversion');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "conversion" VALUES('IGNF','CONV_%s','Conversion for %s',NULL,NULL,'EPSG','1262',NULL,NULL,'PROJ mill',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);""" % (code, escape_literal(d['title']))
+        all_sql.append(sql)
+
+        #sql = """INSERT INTO "crs" VALUES('IGNF','%s','projected');""" % (code)
+        #all_sql.append(sql)
+
+        sql = """INSERT INTO "projected_crs" VALUES('IGNF','%s','%s',NULL,NULL,'EPSG','4499','IGNF','%s','IGNF','CONV_%s','EPSG','1262',NULL,0);""" % (code, escape_literal(d['title']), geogcrs_code, code)
+        all_sql.append(sql)
+
+        all_sql.append('')
+
+    else:
+        assert False, line
+
+
+all_sql.append("""INSERT INTO grid_alternatives(original_grid_name,
+                              proj_grid_name,
+                              proj_grid_format,
+                              proj_method,
+                              inverse_direction,
+                              package_name,
+                              url, direct_download, open_license, directory)
+                      VALUES ('ntf_r93.gsb',    -- as referenced by the IGNF registry
+                              'ntf_r93.gsb',
+                              'NTv2',
+                              'hgridshift',
+                              0,
+                              'proj-datumgrid',
+                              NULL, NULL, NULL, NULL);
+""")
+
+#for k in map_geog_crs:
+#    print k, map_geog_crs[k]
+
+script_dir_name = os.path.dirname(os.path.realpath(__file__))
+sql_dir_name = os.path.join(os.path.dirname(script_dir_name), 'data', 'sql')
+
+f = open(os.path.join(sql_dir_name, 'ignf') + '.sql', 'wb')
+f.write("--- This file has been generated by scripts/build_db_create_ignf.py from the PROJ4 IGNF definition file. DO NOT EDIT !\n\n".encode('UTF-8'))
+for sql in all_sql:
+    f.write((sql + '\n').encode('UTF-8'))
+f.close()
+