Commit f596a8b7 authored by Daniel Scheffler's avatar Daniel Scheffler Committed by Daniel Scheffler
Browse files

Fixed EPSG2WKT returning None in case GDAL_DATA environment variable is not...

Fixed EPSG2WKT returning None in case GDAL_DATA environment variable is not set. Added Test_EPSG2WKT.
parent 40dbf777
Pipeline #1072 passed with stages
in 40 seconds
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import re
import os
import pyproj
......@@ -22,8 +21,7 @@ except ImportError:
from gdalconst import *
def get_proj4info(ds=None,proj=None):
def get_proj4info(ds=None, proj=None):
# type: (gdal.Dataset,str) -> str
"""Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectivly,
e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
......@@ -36,7 +34,7 @@ def get_proj4info(ds=None,proj=None):
if isinstance(proj, str) and proj.startswith('epsg:'):
proj = EPSG2WKT(int(proj.split(':')[1]))
elif isinstance(proj,int):
elif isinstance(proj, int):
proj = EPSG2WKT(proj)
srs.ImportFromWkt(proj)
......@@ -44,12 +42,12 @@ def get_proj4info(ds=None,proj=None):
def proj4_to_dict(proj4):
#type: (str) -> dict
# type: (str) -> dict
"""Converts a PROJ4-like string into a dictionary.
:param proj4: <str> the PROJ4-like string
"""
items = [item for item in proj4.replace('+','').split(' ') if '=' in item]
return {k:v for k,v in [kv.split('=') for kv in items]}
items = [item for item in proj4.replace('+', '').split(' ') if '=' in item]
return {k: v for k, v in [kv.split('=') for kv in items]}
def dict_to_proj4(proj4dict):
......@@ -80,15 +78,17 @@ def prj_equal(prj1, prj2):
if prj1 is None and prj2 is None or prj1 == prj2:
return True
else:
return get_proj4info(proj=prj1)==get_proj4info(proj=prj2)
return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
def isProjectedOrGeographic(prj):
# type: (Union[str, int, dict]) -> Union[str, None]
"""
:param prj: accepts EPSG, Proj4 and WKT projections#
:param prj: accepts EPSG, Proj4 and WKT projections
"""
if prj is None: return None
if prj is None:
return None
srs = osr.SpatialReference()
if prj.startswith('EPSG:'):
srs.ImportFromEPSG(int(prj.split(':')[1]))
......@@ -114,9 +114,21 @@ def EPSG2Proj4(EPSG_code):
def EPSG2WKT(EPSG_code):
# type: (int) -> str
if EPSG_code is not None:
# GDAL_DATA_tmp = None
# if 'GDAL_DATA' in os.environ:
# GDAL_DATA_tmp = os.environ['GDAL_DATA']
# del os.environ['GDAL_DATA']
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG_code)
return srs.ExportToWkt()
wkt = srs.ExportToWkt()
# if GDAL_DATA_tmp:
# os.environ['GDAL_DATA'] = GDAL_DATA_tmp
if not wkt:
raise EnvironmentError(gdal.GetLastErrorMsg())
return wkt
else:
return ''
def _find_epsgfile():
......@@ -139,6 +151,7 @@ def _find_epsgfile():
def WKT2EPSG(wkt, epsgfile=''):
# type: (str) -> Union[int, None]
""" Transform a WKT string to an EPSG code
:param wkt: WKT definition
:param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
......@@ -152,21 +165,21 @@ def WKT2EPSG(wkt, epsgfile=''):
# FIXME PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],
# FIXME UNIT["Meter",1.0]]}
if not isinstance(wkt,str):
raise TypeError("'wkt' must be a string. Received %s." %type(wkt))
code = None #default
if not isinstance(wkt, str):
raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
code = None # default
if not wkt:
return None
p_in = osr.SpatialReference()
s = p_in.ImportFromWkt(wkt)
s = p_in.ImportFromWkt(wkt)
if s == 5: # invalid WKT
raise Exception('Received an invalid WKT string: %s' %wkt)
raise Exception('Received an invalid WKT string: %s' % wkt)
if p_in.IsLocal():
raise Exception('The given WKT is a local coordinate system.')
cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
p_in.AutoIdentifyEPSG()
an = p_in.GetAuthorityName(cstype)
assert an in [None,'EPSG'], "No EPSG code found. Found %s instead." %an
assert an in [None, 'EPSG'], "No EPSG code found. Found %s instead." % an
ac = p_in.GetAuthorityCode(cstype)
if ac is None: # try brute force approach by grokking proj epsg definition file
p_out = p_in.ExportToProj4()
......@@ -179,15 +192,15 @@ def WKT2EPSG(wkt, epsgfile=''):
if m:
code = m.group(1)
break
code = int(code) if code else None # match or no match
code = int(code) if code else None # match or no match
else:
code = int(ac)
return code
def get_UTMzone(ds=None,prj=None):
def get_UTMzone(ds=None, prj=None):
assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
if isProjectedOrGeographic(prj)=='projected':
if isProjectedOrGeographic(prj) == 'projected':
srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjection() if ds else prj)
return srs.GetUTMZone()
......@@ -196,11 +209,11 @@ def get_UTMzone(ds=None,prj=None):
def get_prjLonLat(fmt='wkt'):
# type: (str) -> Any
# type: (str) -> Union[str, dict]
"""Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
:param fmt: <str> target format - 'WKT' or 'PROJ4'
"""
assert re.search('wkt',fmt,re.I) or re.search('Proj4',fmt,re.I), 'unsupported output format'
assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
return srs.ExportToWkt() if re.search('wkt',fmt,re.I) else srs.ExportToProj4()
return srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()
......@@ -11,39 +11,51 @@ Tests for `py_tools_ds.geo.projection` module.
import unittest
from py_tools_ds.geo.projection import WKT2EPSG
from py_tools_ds.geo.projection import WKT2EPSG, EPSG2WKT
wkt_utm = \
"""
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84", 6378137, 298.257223563,
AUTHORITY["EPSG", "7030"]],
AUTHORITY["EPSG", "6326"]],
PRIMEM["Greenwich", 0,
AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.0174532925199433,
AUTHORITY["EPSG", "9122"]],
AUTHORITY["EPSG", "4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin", 0],
PARAMETER["central_meridian", 15],
PARAMETER["scale_factor", 0.9996],
PARAMETER["false_easting", 500000],
PARAMETER["false_northing", 0],
UNIT["metre", 1,
AUTHORITY["EPSG", "9001"]],
AXIS["Easting", EAST],
AXIS["Northing", NORTH],
AUTHORITY["EPSG", "32633"]]
"""
class Test_WKT2EPSG(unittest.TestCase):
def setUp(self):
self.wkt_utm = \
"""
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84", 6378137, 298.257223563,
AUTHORITY["EPSG", "7030"]],
AUTHORITY["EPSG", "6326"]],
PRIMEM["Greenwich", 0,
AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.0174532925199433,
AUTHORITY["EPSG", "9122"]],
AUTHORITY["EPSG", "4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin", 0],
PARAMETER["central_meridian", 15],
PARAMETER["scale_factor", 0.9996],
PARAMETER["false_easting", 500000],
PARAMETER["false_northing", 0],
UNIT["metre", 1,
AUTHORITY["EPSG", "9001"]],
AXIS["Easting", EAST],
AXIS["Northing", NORTH],
AUTHORITY["EPSG", "32633"]]
"""
self.wkt_utm = wkt_utm
def test_UTM_wkt(self):
epsg = WKT2EPSG(self.wkt_utm, epsgfile='')
print('epsg', epsg)
self.assertTrue(isinstance(epsg, int))
class Test_EPSG2WKT(unittest.TestCase):
def setUp(self):
self.epsg_utm = 32636
def test_UTM_epsg(self):
wkt = EPSG2WKT(self.epsg_utm)
self.assertTrue(isinstance(wkt, str), "EPSG2WKT returned a %s object instead of a string!" % type(wkt))
self.assertNotEquals(wkt, "", msg="EPSG2WKT returned an empty WKT string!")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment