Commit b934a84f authored by Daniel Scheffler's avatar Daniel Scheffler

Replaced further code by pyproj>2.1 equivalents. Added test_coord_trafo.py.

Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 384b032d
Pipeline #11898 failed with stage
in 6 minutes and 7 seconds
......@@ -67,7 +67,7 @@ test_py_tools_ds_install:
# - conda config --set channel_priority strict # otherwise gdal or libgdal may be installed from defaults channel
# resolve some requirements with conda
- conda install --yes -q -c conda-forge numpy gdal pyproj>2.0.0 shapely geopandas
- conda install --yes -q -c conda-forge numpy gdal pyproj>2.1.0 shapely geopandas
# run installer
- python setup.py install
......
......@@ -62,7 +62,7 @@ Using conda_, the recommended approach is:
# create virtual environment for py_tools_ds, this is optional
conda create -y -q -c conda-forge --name py_tools_ds python=3
source activate py_tools_ds
conda install -c conda-forge numpy gdal pyproj>2.0.0 shapely scikit-image pandas
conda install -c conda-forge numpy gdal pyproj>2.1.0 shapely scikit-image pandas
Then install py_tools_ds using the pip installer:
......
......@@ -73,16 +73,9 @@ def transform_any_prj(prj_src, prj_tgt, x, y):
:param y:
:return:
"""
from pyproj import __version__ as ver
if ver.startswith('2') and int(ver.split('.')[1]) >= 1:
transformer = pyproj.Transformer.from_crs(get_proj4info(proj=prj_src),
get_proj4info(proj=prj_tgt))
x, y = transformer.transform(x, y)
else:
prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
x, y = pyproj.transform(prj_src, prj_tgt, x, y)
return x, y
transformer = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(prj_src),
pyproj.CRS.from_user_input(prj_tgt))
return transformer.transform(x, y)
def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
......@@ -313,8 +306,8 @@ def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
:param prj_src: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
:param prj_tgt: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
"""
project = partial(pyproj.transform, pyproj.Proj(get_proj4info(proj=prj_src)),
pyproj.Proj(get_proj4info(proj=prj_tgt)))
project = partial(pyproj.transform,
pyproj.CRS.from_user_input(prj_src),
pyproj.CRS.from_user_input(prj_tgt))
return transform(project, shapelyGeometry) # apply projection
......@@ -27,6 +27,7 @@ import warnings
import os
from typing import Union # noqa F401 # flake8 issue
from pyproj import CRS
import numpy as np
try:
from osgeo import osr
......@@ -37,7 +38,7 @@ except ImportError:
import gdal
from .coord_trafo import transform_any_prj
from .projection import get_proj4info, proj4_to_dict, isLocal
from .projection import isLocal
__author__ = "Daniel Scheffler"
......@@ -130,14 +131,17 @@ class Geocoding(object):
if isLocal(prj):
self.prj_name = 'Arbitrary'
else:
# get prj_name and datum
proj4 = proj4_to_dict(get_proj4info(proj=prj)) # type: dict
self.prj_name = \
'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \
'UTM' if proj4['proj'] == 'utm' else proj4['proj']
# FIXME there is no 'datum' key in case of, e.g., LAEA projection
self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum'] # proj4['ellps']?
self.units = proj4['unit'] if 'unit' in proj4 else self.units
with warnings.catch_warnings():
warnings.filterwarnings("ignore", message="You will likely lose important projection information")
# get prj_name and datum
proj4 = CRS(prj).to_dict() # FIXME avoid to convert to proj4
self.prj_name = \
'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \
'UTM' if proj4['proj'] == 'utm' else proj4['proj']
# FIXME there is no 'datum' key in case of, e.g., LAEA projection # still true?
self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum'] # proj4['ellps']?
self.units = proj4['unit'] if 'unit' in proj4 else self.units
if self.prj_name == 'UTM':
self.utm_zone = srs.GetUTMZone()
......
......@@ -42,25 +42,15 @@ __author__ = "Daniel Scheffler"
gdal_env.try2set_GDAL_DATA()
def get_proj4info(ds=None, proj=None):
# type: (gdal.Dataset, Union[str, int]) -> str
"""Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectively,
def get_proj4info(proj=None):
# type: (Union[str, int]) -> str
"""Returns PROJ4 formatted projection info for the given projection.
e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
:param ds: <gdal.Dataset> the gdal dataset to get PROJ4 info for
:param proj: <str,int> the projection to get PROJ4 formatted info for (WKT or 'epsg:1234' or <EPSG_int>)
"""
assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
srs = osr.SpatialReference()
proj = ds.GetProjection() if ds else proj
if isinstance(proj, str) and proj.startswith('epsg:'):
proj = EPSG2WKT(int(proj.split(':')[1]))
elif isinstance(proj, int):
proj = EPSG2WKT(proj)
srs.ImportFromWkt(proj)
return srs.ExportToProj4().strip()
return CRS.from_user_input(proj).to_proj4().strip()
def proj4_to_dict(proj4):
......@@ -97,7 +87,10 @@ 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)
crs1 = CRS.from_user_input(prj1)
crs2 = CRS.from_user_input(prj2)
return crs1.equals(crs2)
def isProjectedOrGeographic(prj):
......
......@@ -36,7 +36,7 @@ version = {}
with open("py_tools_ds/version.py") as version_file:
exec(version_file.read(), version)
requirements = ['gdal>=2.1.0', 'numpy', 'shapely', 'six', 'pandas', 'scikit-image>=0.16.0', 'geopandas', 'pyproj>2.0.0',
requirements = ['gdal>=2.1.0', 'numpy', 'shapely', 'six', 'pandas', 'scikit-image>=0.16.0', 'geopandas', 'pyproj>2.1.0',
'spectral']
setup_requirements = [] # TODO(danschef): put setup requirements (distutils extensions, etc.) here
test_requirements = requirements + ["coverage", "nose", "nose2", "nose-htmloutput", "rednose"]
......
......@@ -10,7 +10,7 @@ dependencies:
- gdal>=2.1.0
- conda-forge::libgdal # force to use conda-forge for libgdal to avoid package version incompatiblies due to mixed channels (libkea.so.1.4.7: cannot open shared object file: No such file or directory)
# - kealib # fix for libkea.so.1.4.7: cannot open shared object file: No such file or directory (not needed as long as libgdal and gdal are conda-forge packages)
- pyproj>2.0.0
- pyproj>2.1.0
- scikit-image>=0.16.0
- geopandas
- pip:
......
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# py_tools_ds
#
# Copyright (C) 2019 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
"""
test_coord_trafo
----------------
Tests for `py_tools_ds.geo.coord_trafo` module.
"""
import unittest
from shapely.geometry import Polygon
from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry
poly_local = Polygon([(5708.2, -3006), (5708, -3262), (5452, -3262), (5452, -3006), (5708, -3006)])
class Test_reproject_shapelyGeometry(unittest.TestCase):
def test_reproject_shapelyGeometry(self):
poly_lonlat = reproject_shapelyGeometry(poly_local, 32636, 4326)
self.assertTrue(isinstance(poly_lonlat, Polygon))
......@@ -71,6 +71,7 @@ wkt_utm = \
AXIS["Northing", NORTH],
AUTHORITY["EPSG", "32633"]]
"""
wkt_utm = ' '.join(wkt_utm.split())
class Test_geotransform2mapinfo(unittest.TestCase):
......
......@@ -65,7 +65,7 @@ wkt_utm = \
class Test_WKT2EPSG(unittest.TestCase):
def setUp(self):
self.wkt_utm = wkt_utm
self.wkt_utm = ' '.join(wkt_utm.split())
def test_UTM_wkt(self):
epsg = WKT2EPSG(self.wkt_utm)
......
......@@ -64,6 +64,7 @@ class Test_boxObj(TestCase):
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32636"]]
"""
self.prj = ' '.join(self.prj.split())
self.box = boxObj(wp=self.wp,
ws=self.ws,
......@@ -175,7 +176,7 @@ class Test_boxObj(TestCase):
self.box.get_coordArray_MapXY()
# test to return coords array in another projection
self.box.get_coordArray_MapXY(prj="""
self.box.get_coordArray_MapXY(prj=' '.join("""
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
......@@ -198,7 +199,8 @@ class Test_boxObj(TestCase):
AXIS["Easting", EAST],
AXIS["Northing", NORTH],
AUTHORITY["EPSG", "32633"]]
""")
""".split())
)
box = self.box
box.prj = None
......
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