# BEGIN OF LICENSE NOTE
# This file is part of Pyoints.
# Copyright (c) 2018, Sebastian Lamprecht, Trier University,
# lamprecht@uni-trier.de
#
# Pyoints is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Pyoints is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Pyoints. If not, see <https://www.gnu.org/licenses/>.
# END OF LICENSE NOTE
"""Coordinate Reference Systems and two dimensional geograpic coordinate
transformations.
"""
import pyproj
import numpy as np
from osgeo import osr
from . import (
assertion,
)
from .misc import print_rounded
# Global proj4 definitions
WGS84 = '+proj=latlong +datum=WGS84 +to +proj=latlong +datum=WGS84 +units=m ' \
'+no_defs'
[docs]class Proj:
"""Wrapper class for different commonly coordinate reference system
formats.
Parameters
----------
proj4 : optional, str
Coordinate reference system definition in Proj4 format. WGS84, if None
or empty string.
Attributes
----------
proj4 : str
Projection in Proj4 format.
wkt : str
Projection in Well Known Text format.
pyproj : `pyproj.Proj`
Projection as `pyproj.Proj` instance.
Examples
--------
Create from EPSG code.
>>> proj = Proj.from_epsg(4326)
>>> print('AUTHORITY["EPSG","4326"]' in proj.wkt)
True
>>> print('WGS84' in proj.proj4)
True
Create from Proj4 string.
>>> proj = Proj.from_proj4('+proj=longlat +datum=WGS84 +no_defs')
>>> print(proj.proj4)
+proj=longlat +datum=WGS84 +no_defs
>>> print('AUTHORITY["EPSG","4326"]' in proj.wkt)
True
Create from Well Known Text.
>>> proj = Proj.from_wkt(proj.wkt)
>>> print('WGS84' in proj.proj4)
True
"""
def __init__(self, proj4=WGS84):
if proj4 is None or not isinstance(proj4, str) or proj4 == '':
raise ValueError("proj4 not defined")
self._proj4 = proj4.rstrip()
@property
def proj4(self):
return self._proj4
@property
def wkt(self):
sr = osr.SpatialReference()
sr.ImportFromProj4(self.proj4)
return sr.ExportToWkt()
@property
def pyproj(self):
return pyproj.Proj(self.proj4)
@property
def osr(self):
srs = osr.SpatialReference()
srs.ImportFromProj4(self.proj4)
return srs
def __str__(self):
return 'proj4: %s' % str(self.proj4)
[docs] @classmethod
def from_proj4(cls, proj4):
"""Create `Proj` object from Proj4 format.
Parameters
----------
proj4 : str
Coordinate projection definition in Proj4 format.
"""
return Proj(proj4)
[docs] @classmethod
def from_wkt(cls, wkt):
"""Create `Proj` object from Well Known Text.
Parameters
----------
wkt : str
Coordinate projection definition in Well Known Text format.
"""
if not isinstance(wkt, str):
raise TypeError("'wkt' needs to be a string")
proj4 = osr.SpatialReference(wkt=wkt).ExportToProj4()
if proj4 == '':
raise ValueError("WKT unknown")
return Proj.from_proj4(proj4)
[docs] @classmethod
def from_epsg(cls, epsg):
"""Create `Proj` object from EPSG code.
Parameters
----------
epsg : int
Coordinate projection definition in EPSG format.
"""
if not isinstance(epsg, int):
raise TypeError("'epsg' needs to be an integer")
sr = osr.SpatialReference()
sr.ImportFromEPSG(epsg)
proj4 = sr.ExportToProj4()
if proj4 == '':
raise ValueError("epsg code '%i' unknown" % epsg)
return Proj.from_proj4(proj4)
[docs]def project(coords, from_proj, to_proj):
"""Applies the coordinate transformation.
Parameters
----------
coords : array_like(Number, shape=(n, k))
Represents `n` points of `k` dimensions to transform.
from_proj,to_proj : `Proj`
Define the coordinate transformation from the origin projection system
`from_proj` to the target projection system `to_proj`.
See Also
--------
GeoTransform
Returns
-------
np.ndarray(Number, shape=(n, k))
Transformed coordinates.
"""
geoTransform = GeoTransform(from_proj, to_proj)
return geoTransform(coords)
