JezK

Edit File: epsg_tr.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ******************************************************************************
#  $Id$
#
#  Project:  CFS OGC MapServer
#  Purpose:  Script to create WKT and PROJ.4 dictionaries for EPSG GCS/PCS
#            codes.
#  Author:   Frank Warmerdam, warmerdam@pobox.com
#
# ******************************************************************************
#  Copyright (c) 2001, Frank Warmerdam
#  Copyright (c) 2009-2010, 2019, Even Rouault <even dot rouault at spatialys.com>
#  Copyright (c) 2021, Idan Miara <idan@miara.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 sys
from typing import Optional

from osgeo import gdal, osr
from osgeo_utils.auxiliary.gdal_argparse import GDALArgumentParser, GDALScript

def Usage():
    print(
        f"Usage: {sys.argv[0]} -- This is a sample. Read source to know how to use. --"
    )
    return 2

def trHandleCode(set_srid, srs, auth_name, code, deprecated, output_format):
    if output_format == "-pretty_wkt":
        print("%s:%s" % (auth_name, str(code)))
        print(srs.ExportToPrettyWkt())

if output_format == "-xml":
        print(srs.ExportToXML())

if output_format == "-wkt":
        print(f"EPSG:{code}")
        print(srs.ExportToWkt())

if output_format == "-proj4":
        out_string = srs.ExportToProj4()

name = srs.GetName()

print("# %s" % name)
        if out_string.find("+proj=") > -1:
            print("<%s> %s <>" % (str(code), out_string))
        else:
            print(
                "# Unable to translate coordinate system "
                "%s:%s into PROJ.4 format." % (auth_name, str(code))
            )
            print("#")

if output_format == "-postgis":

if code in set_srid:
            if auth_name == "ESRI":
                if int(code) < 32767:
                    return
        assert code not in set_srid, (auth_name, code)
        set_srid.add(code)

name = srs.GetName()
        if deprecated and "deprecated" not in name:
            name += " (deprecated)"
        wkt = srs.ExportToWkt()
        proj4text = srs.ExportToProj4()

print("---")
        print("--- %s %s : %s" % (auth_name, str(code), name))
        print("---")

if proj4text is None or len(proj4text) == 0:
            print("-- (unable to translate to PROJ.4)")
        else:
            wkt = gdal.EscapeString(wkt, scheme=gdal.CPLES_SQL)
            proj4text = gdal.EscapeString(proj4text, scheme=gdal.CPLES_SQL)
            print(
                'INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (%d,\'%s\',%d,\'%s\',\'%s\');'
                % (int(code), auth_name, int(code), wkt, proj4text)
            )

# INGRES COPY command input.
    if output_format == "-copy":

try:
            wkt = srs.ExportToWkt()
            proj4text = srs.ExportToProj4()

print(
                "%s\t%d%s\t%s\t%d%s\t%d%s\n"
                % (
                    str(code),
                    4,
                    auth_name,
                    str(code),
                    len(wkt),
                    wkt,
                    len(proj4text),
                    proj4text,
                )
            )
        except Exception:
            pass

def epsg_tr(output_format: str = "-pretty_wkt", authority: Optional[str] = None):
    # Output BEGIN transaction for PostGIS
    if output_format == "-postgis":
        print("BEGIN;")

# loop over all codes to generate output

if authority:
        authorities = [authority]
    elif output_format == "-postgis":
        authorities = ["EPSG", "ESRI"]
    else:
        authorities = ["EPSG", "ESRI", "IGNF"]

set_srid = set()
    for authority in authorities:
        if authority in ("EPSG", "ESRI"):
            set_codes_geographic = set()
            set_codes_geographic_3d = set()
            set_codes_projected = set()
            set_codes_geocentric = set()
            set_codes_compound = set()
            set_deprecated = set()

for crs_info in osr.GetCRSInfoListFromDatabase(authority):
                code = int(crs_info.code)
                if crs_info.type == osr.OSR_CRS_TYPE_COMPOUND:
                    set_codes_compound.add(code)
                elif crs_info.type == osr.OSR_CRS_TYPE_GEOGRAPHIC_3D:
                    set_codes_geographic_3d.add(code)
                elif crs_info.type == osr.OSR_CRS_TYPE_GEOGRAPHIC_2D:
                    set_codes_geographic.add(code)
                elif crs_info.type == osr.OSR_CRS_TYPE_PROJECTED:
                    set_codes_projected.add(code)
                elif crs_info.type == osr.OSR_CRS_TYPE_GEOCENTRIC:
                    set_codes_geocentric.add(code)

if crs_info.deprecated:
                    set_deprecated.add(code)

set_codes_geographic = sorted(set_codes_geographic)
            set_codes_geographic_3d = sorted(set_codes_geographic_3d)
            set_codes_projected = sorted(set_codes_projected)
            set_codes_geocentric = sorted(set_codes_geocentric)
            set_codes_compound = sorted(set_codes_compound)
            for typestr, set_codes in (
                ("Geographic 2D CRS", set_codes_geographic),
                ("Projected CRS", set_codes_projected),
                ("Geocentric CRS", set_codes_geocentric),
                ("Compound CRS", set_codes_compound),
                ("Geographic 3D CRS", set_codes_geographic_3d),
            ):
                if set_codes and output_format == "-postgis":
                    print("-" * 80)
                    print("--- " + authority + " " + typestr)
                    print("-" * 80)

for code in set_codes:
                    srs = osr.SpatialReference()
                    srs.SetFromUserInput(authority + ":" + str(code))
                    deprecated = False
                    if code in set_deprecated:
                        deprecated = True
                    trHandleCode(
                        set_srid, srs, authority, str(code), deprecated, output_format
                    )

else:
            for crs_info in osr.GetCRSInfoListFromDatabase(authority):
                srs = osr.SpatialReference()
                srs.SetFromUserInput(authority + ":" + crs_info.code)
                trHandleCode(
                    set_srid,
                    srs,
                    authority,
                    crs_info.code,
                    crs_info.deprecated,
                    output_format,
                )

# Output COMMIT transaction for PostGIS
    if output_format == "-postgis":
        print("COMMIT;")
        print("VACUUM ANALYZE spatial_ref_sys;")
    return 0

class EPSG_Table(GDALScript):
    def __init__(self):
        super().__init__()
        self.title = "Create WKT and PROJ.4 dictionaries for EPSG GCS/PCS codes"

def get_parser(self, argv) -> GDALArgumentParser:
        parser = self.parser

parser.add_argument(
            "-authority",
            dest="authority",
            metavar="name",
            type=str,
            help="Authority name",
        )

group = parser.add_mutually_exclusive_group()
        group.add_argument(
            "-wkt",
            const="wkt",
            dest="output_format",
            action="store_const",
            help="Well Known Text",
        )
        group.add_argument(
            "-pretty_wkt",
            const="pretty_wkt",
            dest="output_format",
            action="store_const",
            help="Pretty Well Known Text",
        )
        group.add_argument(
            "-proj4",
            const="proj4",
            dest="output_format",
            action="store_const",
            help="proj string",
        )
        group.add_argument(
            "-postgis",
            const="postgis",
            dest="output_format",
            action="store_const",
            help="postgis",
        )
        group.add_argument(
            "-xml", const="xml", dest="output_format", action="store_const", help="XML"
        )
        group.add_argument(
            "-copy",
            const="copy",
            dest="output_format",
            action="store_const",
            help="Table",
        )

return parser

def augment_kwargs(self, kwargs) -> dict:
        if not kwargs.get("output_format"):
            kwargs["output_format"] = "pretty_wkt"
        kwargs["output_format"] = "-" + kwargs["output_format"]
        return kwargs

def doit(self, **kwargs):
        return epsg_tr(**kwargs)

def main(argv=sys.argv):
    if len(sys.argv) < 2:
        return Usage()
    return EPSG_Table().main(argv)

if __name__ == "__main__":
    sys.exit(main(sys.argv))

JezK

Edit File: epsg_tr.py

FAIRNETZT