Database: import common projections from ESRI projected CRS in structured form

That is Transverse_Mercator/Gauss_Kruger, Lambert_Conformal_Conic and
Hotine_Oblique_Mercator_Azimuth_Natural_Origin

Decreases proj.db from 5 853 184 bytes to 5 189 632 bytes.

1 files changed, 264 insertions, 4 deletions

diff --git a/scripts/build_db_from_esri.py b/scripts/build_db_from_esri.py
index af0a67b5..6de2eb85 100755
--- a/scripts/build_db_from_esri.py
+++ b/scripts/build_db_from_esri.py
@@ -615,9 +615,122 @@ def import_geogcs():
 
 ########################
 
+def parse_wkt(s, level):
+    if s[0] == '"':
+        return s
+    pos = s.find('[')
+    if pos < 0:
+        return s
+    return { s[0:pos] : parse_wkt_array(s[pos+1:-1], level + 1) }
+
+def parse_wkt_array(s, level = 0):
+    ar = []
+    in_string = False
+    cur_token = ''
+    indent_level = 0
+    for c in s:
+        if in_string:
+            if c == '"':
+                in_string = False
+            cur_token += c
+        elif c == '"':
+            cur_token += c
+            in_string = True
+        elif c == '[':
+            cur_token += c
+            indent_level += 1
+        elif c == ']':
+            cur_token += c
+            indent_level -= 1
+            assert indent_level >= 0
+        elif indent_level == 0 and  c == ',':
+            ar.append(parse_wkt(cur_token, level + 1))
+            cur_token = ''
+        else:
+            cur_token += c
+    assert indent_level == 0
+    if cur_token:
+        ar.append(parse_wkt(cur_token, level + 1))
+
+    if level == 0:
+        d = {}
+        for elt in ar:
+            assert type(elt) == type({})
+            assert len(elt) == 1
+            if 'PROJECTION' in elt:
+                assert len(elt['PROJECTION']) == 1, elt['PROJECTION']
+                assert 'PROJECTION' not in d
+                name = elt['PROJECTION'][0]
+                assert name[0] == '"' and name[-1] == '"', name
+                name = name[1:-1]
+                d['PROJECTION'] = name
+            elif 'PARAMETER' in elt:
+                assert len(elt['PARAMETER']) == 2, elt['PARAMETER']
+                name = elt['PARAMETER'][0]
+                assert name[0] == '"' and name[-1] == '"', name
+                name = name[1:-1]
+                assert name not in d
+                d[name] = elt['PARAMETER'][1]
+            elif 'UNIT' in elt:
+                assert len(elt['UNIT']) == 2, elt['UNIT']
+                name = elt['UNIT'][0]
+                assert name[0] == '"' and name[-1] == '"', name
+                name = name[1:-1]
+                assert 'UNIT_NAME' not in d
+                d['UNIT_NAME'] = name
+                d['UNIT_VALUE'] = elt['UNIT'][1]
+            else:
+                assert True
+        return d
+    else:
+        return ar
 
-map_projcs_esri_name_to_auth_code = {}
+########################
+
+def get_cs(parsed_conv_wkt):
+
+    UNIT_NAME = parsed_conv_wkt['UNIT_NAME']
+    UNIT_VALUE = parsed_conv_wkt['UNIT_VALUE']
+
+    if UNIT_NAME == 'Meter':
+        uom_code = '9001'
+        cs_auth_name = 'EPSG'
+        cs_code = '4400'
+        assert UNIT_VALUE == '1.0', UNIT_VALUE
+    elif UNIT_NAME == 'Foot':
+        uom_code = '9002'
+        cs_auth_name = 'ESRI'
+        cs_code = UNIT_NAME
+        assert UNIT_VALUE == '0.3048', UNIT_VALUE
+    elif UNIT_NAME == 'Foot_US':
+        uom_code = '9003'
+        cs_auth_name = 'ESRI'
+        cs_code = UNIT_NAME
+        assert UNIT_VALUE == '0.3048006096012192', UNIT_VALUE
+    elif UNIT_NAME == 'Yard_Indian_1937':
+        uom_code = '9085'
+        cs_auth_name = 'ESRI'
+        cs_code = UNIT_NAME
+        assert UNIT_VALUE == '0.91439523', UNIT_VALUE
+    else:
+        assert False, UNIT_NAME
+
+    if cs_auth_name == 'ESRI' and cs_code not in set_esri_cs_code:
+        sql = """INSERT INTO "coordinate_system" VALUES('ESRI','%s','Cartesian',2);""" % cs_code
+        all_sql.append(sql)
+        sql = """INSERT INTO "axis" VALUES('ESRI','%d','Easting','E','east','ESRI','%s',1,'EPSG','%s');""" % (2 * len(set_esri_cs_code) + 1, cs_code, uom_code)
+        all_sql.append(sql)
+        sql = """INSERT INTO "axis" VALUES('ESRI','%d','Northing','N','north','ESRI','%s',2,'EPSG','%s');""" % (2 * len(set_esri_cs_code) + 2, cs_code, uom_code)
+        all_sql.append(sql)
+        set_esri_cs_code.add(cs_code)
+
+    return cs_auth_name, cs_code, uom_code
+
+########################
 
+map_projcs_esri_name_to_auth_code = {}
+set_esri_cs_code = set()
+map_conversion_sql_to_code = {}
 
 def import_projcs():
     with open(os.path.join(path_to_csv, 'pe_list_projcs.csv'), 'rt') as csvfile:
@@ -705,6 +818,10 @@ def import_projcs():
                 end_pos += pos
                 geogcs_name = wkt[pos:end_pos]
 
+                pos = wkt.find('PROJECTION[')
+                assert pos >= 0
+                parsed_conv_wkt = parse_wkt_array(wkt[pos:-1])
+
                 assert geogcs_name in map_geogcs_esri_name_to_auth_code, (
                     geogcs_name, row)
                 geogcs_auth_name, geogcs_code = map_geogcs_esri_name_to_auth_code[geogcs_name]
@@ -716,9 +833,152 @@ def import_projcs():
 
                 deprecated = 1 if row[idx_deprecated] == 'yes' else 0
 
-                sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,NULL,NULL,'%s','%s',NULL,NULL,'%s','%s','%s',%d);""" % (
-                    code, esri_name, geogcs_auth_name, geogcs_code, area_auth_name, area_code, escape_literal(wkt), deprecated)
-                all_sql.append(sql)
+                method = parsed_conv_wkt['PROJECTION']
+
+                if 'UNIT["Degree",' in wkt:
+                    ang_uom_code = '9102'
+                elif 'UNIT["Grad",' in wkt:
+                    ang_uom_code = '9105'
+                else:
+                    assert False, wkt
+
+                if method in ('Transverse_Mercator', 'Gauss_Kruger'):
+                    assert len(parsed_conv_wkt) == 1 + 5 + 2
+                    False_Easting = parsed_conv_wkt['False_Easting']
+                    False_Northing = parsed_conv_wkt['False_Northing']
+                    Central_Meridian = parsed_conv_wkt['Central_Meridian']
+                    Scale_Factor = parsed_conv_wkt['Scale_Factor']
+                    Latitude_Of_Origin = parsed_conv_wkt['Latitude_Of_Origin']
+
+                    cs_auth_name, cs_code, uom_code = get_cs(parsed_conv_wkt)
+
+                    conv_name = 'unnamed'
+                    if method == 'Gauss_Kruger' and 'GK_' not in esri_name and 'Gauss' not in esri_name:
+                        conv_name = esri_name + " (Gauss Kruger)"
+
+                    cursor.execute(
+                        """SELECT code FROM conversion WHERE auth_name = 'EPSG' AND
+                        method_code = '9807' AND
+                        param1_code = '8801' AND param1_value = ? AND param1_uom_code = ? AND
+                        param2_code = '8802' AND param2_value = ? AND param2_uom_code = ? AND
+                        param3_code = '8805' AND param3_value = ? AND param3_uom_code = '9201' AND
+                        param4_code = '8806' AND param4_value = ? AND param4_uom_code = ? AND
+                        param5_code = '8807' AND param5_value = ? AND param5_uom_code = ?""", (Latitude_Of_Origin, ang_uom_code, Central_Meridian, ang_uom_code, Scale_Factor, False_Easting, uom_code, False_Northing, uom_code))
+                    src_row = cursor.fetchone()
+                    if conv_name == 'unnamed' and src_row:
+                        conv_auth_name = 'EPSG'
+                        conv_code = src_row[0]
+                    else:
+                        conv_auth_name = 'ESRI'
+                        conv_code = code
+
+                        sql = """INSERT INTO "conversion" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','EPSG','9807','Transverse Mercator','EPSG','8801','Latitude of natural origin',%s,'EPSG','%s','EPSG','8802','Longitude of natural origin',%s,'EPSG','%s','EPSG','8805','Scale factor at natural origin',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','%s','EPSG','8807','False northing',%s,'EPSG','%s',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,%d);""" % (
+                            code, conv_name, area_auth_name, area_code, Latitude_Of_Origin, ang_uom_code, Central_Meridian, ang_uom_code, Scale_Factor, False_Easting, uom_code, False_Northing, uom_code, deprecated)
+
+                        sql_extract = sql[sql.find('NULL,NULL'):]
+                        if conv_name != 'unnamed' or sql_extract not in map_conversion_sql_to_code:
+                            all_sql.append(sql)
+                            map_conversion_sql_to_code[sql_extract] = conv_code
+                        else:
+                            conv_code = map_conversion_sql_to_code[sql_extract]
+
+                    sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','%s','%s','%s','%s','%s','%s',NULL,%d);""" % (
+                        code, esri_name, cs_auth_name, cs_code, geogcs_auth_name, geogcs_code, conv_auth_name, conv_code, area_auth_name, area_code, deprecated)
+                    all_sql.append(sql)
+
+                elif method == 'Hotine_Oblique_Mercator_Azimuth_Natural_Origin':
+                    assert len(parsed_conv_wkt) == 1 + 6 + 2
+                    False_Easting = parsed_conv_wkt['False_Easting']
+                    False_Northing = parsed_conv_wkt['False_Northing']
+                    Scale_Factor = parsed_conv_wkt['Scale_Factor']
+                    Azimuth = parsed_conv_wkt['Azimuth']
+                    Longitude_Of_Center = parsed_conv_wkt['Longitude_Of_Center']
+                    Latitude_Of_Center = parsed_conv_wkt['Latitude_Of_Center']
+
+                    cs_auth_name, cs_code, uom_code = get_cs(parsed_conv_wkt)
+
+                    conv_name = 'unnamed'
+                    conv_auth_name = 'ESRI'
+                    conv_code = code
+
+                    sql = """INSERT INTO "conversion" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','EPSG','9812','Hotine Oblique Mercator (variant A)','EPSG','8811','Latitude of projection centre',%s,'EPSG','%s','EPSG','8812','Longitude of projection centre',%s,'EPSG','%s','EPSG','8813','Azimuth of initial line',%s,'EPSG','%s','EPSG','8814','Angle from Rectified to Skew Grid',%s,'EPSG','%s','EPSG','8815','Scale factor on initial line',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','%s','EPSG','8807','False northing',%s,'EPSG','%s',%d);""" % (
+                        code, conv_name, area_auth_name, area_code, Latitude_Of_Center, ang_uom_code, Longitude_Of_Center, ang_uom_code, Azimuth, ang_uom_code, Azimuth, ang_uom_code, Scale_Factor, False_Easting, uom_code, False_Northing, uom_code, deprecated)
+
+                    sql_extract = sql[sql.find('NULL,NULL'):]
+                    if conv_name != 'unnamed' or sql_extract not in map_conversion_sql_to_code:
+                        all_sql.append(sql)
+                        map_conversion_sql_to_code[sql_extract] = conv_code
+                    else:
+                        conv_code = map_conversion_sql_to_code[sql_extract]
+
+                    sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','%s','%s','%s','%s','%s','%s',NULL,%d);""" % (
+                        code, esri_name, cs_auth_name, cs_code, geogcs_auth_name, geogcs_code, conv_auth_name, conv_code, area_auth_name, area_code, deprecated)
+                    all_sql.append(sql)
+
+                elif method == 'Lambert_Conformal_Conic' and 'Standard_Parallel_2' in parsed_conv_wkt:
+                    assert len(parsed_conv_wkt) == 1 + 6 + 2
+                    False_Easting = parsed_conv_wkt['False_Easting']
+                    False_Northing = parsed_conv_wkt['False_Northing']
+                    Central_Meridian = parsed_conv_wkt['Central_Meridian']
+                    Standard_Parallel_1 = parsed_conv_wkt['Standard_Parallel_1']
+                    Standard_Parallel_2 = parsed_conv_wkt['Standard_Parallel_2']
+                    Latitude_Of_Origin = parsed_conv_wkt['Latitude_Of_Origin']
+
+                    cs_auth_name, cs_code, uom_code = get_cs(parsed_conv_wkt)
+
+                    conv_name = 'unnamed'
+                    conv_auth_name = 'ESRI'
+                    conv_code = code
+
+                    sql = """INSERT INTO "conversion" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','EPSG','9802','Lambert Conic Conformal (2SP)','EPSG','8821','Latitude of false origin',%s,'EPSG','%s','EPSG','8822','Longitude of false origin',%s,'EPSG','%s','EPSG','8823','Latitude of 1st standard parallel',%s,'EPSG','%s','EPSG','8824','Latitude of 2nd standard parallel',%s,'EPSG','%s','EPSG','8826','Easting at false origin',%s,'EPSG','%s','EPSG','8827','Northing at false origin',%s,'EPSG','%s',NULL,NULL,NULL,NULL,NULL,NULL,%d);""" % (
+                        code, conv_name, area_auth_name, area_code, Latitude_Of_Origin, ang_uom_code, Central_Meridian, ang_uom_code, Standard_Parallel_1, ang_uom_code, Standard_Parallel_2, ang_uom_code, False_Easting, uom_code, False_Northing, uom_code, deprecated)
+
+                    sql_extract = sql[sql.find('NULL,NULL'):]
+                    if conv_name != 'unnamed' or sql_extract not in map_conversion_sql_to_code:
+                        all_sql.append(sql)
+                        map_conversion_sql_to_code[sql_extract] = conv_code
+                    else:
+                        conv_code = map_conversion_sql_to_code[sql_extract]
+
+                    sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','%s','%s','%s','%s','%s','%s',NULL,%d);""" % (
+                        code, esri_name, cs_auth_name, cs_code, geogcs_auth_name, geogcs_code, conv_auth_name, conv_code, area_auth_name, area_code, deprecated)
+                    all_sql.append(sql)
+
+                elif method == 'Lambert_Conformal_Conic' and 'Scale_Factor' in parsed_conv_wkt:
+                    assert len(parsed_conv_wkt) == 1 + 6 + 2
+                    False_Easting = parsed_conv_wkt['False_Easting']
+                    False_Northing = parsed_conv_wkt['False_Northing']
+                    Central_Meridian = parsed_conv_wkt['Central_Meridian']
+                    Standard_Parallel_1 = parsed_conv_wkt['Standard_Parallel_1']
+                    Scale_Factor = parsed_conv_wkt['Scale_Factor']
+                    Latitude_Of_Origin = parsed_conv_wkt['Latitude_Of_Origin']
+                    assert Standard_Parallel_1 == Latitude_Of_Origin
+
+                    cs_auth_name, cs_code, uom_code = get_cs(parsed_conv_wkt)
+
+                    conv_name = 'unnamed'
+                    conv_auth_name = 'ESRI'
+                    conv_code = code
+
+                    sql = """INSERT INTO "conversion" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','EPSG','9801','Lambert Conic Conformal (1SP)','EPSG','8801','Latitude of natural origin',%s,'EPSG','%s','EPSG','8802','Longitude of natural origin',%s,'EPSG','%s','EPSG','8805','Scale factor at natural origin',%s,'EPSG','9201','EPSG','8806','False easting',%s,'EPSG','%s','EPSG','8807','False northing',%s,'EPSG','%s',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,%d);""" % (
+                        code, conv_name, area_auth_name, area_code, Latitude_Of_Origin, ang_uom_code, Central_Meridian, ang_uom_code, Scale_Factor, False_Easting, uom_code, False_Northing, uom_code, deprecated)
+
+                    sql_extract = sql[sql.find('NULL,NULL'):]
+                    if conv_name != 'unnamed' or sql_extract not in map_conversion_sql_to_code:
+                        all_sql.append(sql)
+                        map_conversion_sql_to_code[sql_extract] = conv_code
+                    else:
+                        conv_code = map_conversion_sql_to_code[sql_extract]
+
+                    sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,'%s','%s','%s','%s','%s','%s','%s','%s',NULL,%d);""" % (
+                        code, esri_name, cs_auth_name, cs_code, geogcs_auth_name, geogcs_code, conv_auth_name, conv_code, area_auth_name, area_code, deprecated)
+                    all_sql.append(sql)
+
+                else:
+
+                    sql = """INSERT INTO "projected_crs" VALUES('ESRI','%s','%s',NULL,NULL,NULL,NULL,'%s','%s',NULL,NULL,'%s','%s','%s',%d);""" % (
+                        code, esri_name, geogcs_auth_name, geogcs_code, area_auth_name, area_code, escape_literal(wkt), deprecated)
+                    all_sql.append(sql)
 
                 if deprecated and code != latestWkid:
                     mapDeprecatedToNonDeprecated[code] = latestWkid