Commit fa0b0006 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'bugfix/fix_ETRS_projection_incompatibility' into 'master'

Bugfix/fix etrs projection incompatibility

Closes #7

See merge request !1
parents 5bd37c1e 56610a18
Pipeline #3729 failed with stages
in 16 minutes and 6 seconds
......@@ -2,16 +2,22 @@
import re
import math
import warnings
from tempfile import NamedTemporaryFile
import os
from typing import Union # noqa F401 # flake8 issue
import numpy as np
try:
from osgeo import osr
from osgeo import gdal
except ImportError:
# noinspection PyPackageRequirements
import osr
import gdal
from .coord_trafo import transform_any_prj
from .projection import get_proj4info, proj4_to_dict, isLocal
from ..io.raster.gdal import get_GDAL_ds_inmem
__author__ = "Daniel Scheffler"
......@@ -25,6 +31,7 @@ class Geocoding(object):
:param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
:param prj: GDAL Projection - WKT Format
"""
# FIXME this class only supports UTM and Geographic Lat/Lon projections
self.prj_name = 'Arbitrary'
self.ul_x_map = 0.
self.ul_y_map = 0.
......@@ -49,7 +56,7 @@ class Geocoding(object):
self.from_geotransform_projection(gt, prj)
def from_geotransform_projection(self, gt, prj):
# type: (Union[list, tuple], str) -> self
# type: (Union[list, tuple], str) -> 'Geocoding'
"""Create Geocoding object from GDAL GeoTransform + WKT projection string.
HOW COMPUTATION OF RADIANTS WORKS:
......@@ -108,6 +115,7 @@ class Geocoding(object):
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
......@@ -121,13 +129,12 @@ class Geocoding(object):
return self
def from_mapinfo(self, mapinfo):
# type: (Union[list, tuple]) -> self
# type: (Union[list, tuple]) -> 'Geocoding'
"""Create Geocoding object from ENVI map info.
:param mapinfo: ENVI map info, e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
:return: instance of Geocoding
"""
# type: (Union[list, tuple]) -> self
if mapinfo:
# validate input map info
if not isinstance(mapinfo, (list, tuple)):
......@@ -163,7 +170,7 @@ class Geocoding(object):
return self
def to_geotransform(self):
# type: () -> list
# type: () -> tuple
"""Return GDAL GeoTransform list using the attributes of the Geocoding instance.
For equations, see:
......@@ -181,7 +188,7 @@ class Geocoding(object):
gsd_y, rot_2 = (self.gsd_y, 0) if self.rot_2_deg == 0 else \
(math.cos(self.rot_2_rad) * self.gsd_y, math.sin(self.rot_2_rad) * self.gsd_y)
return [ul_map_x, gsd_x, rot_1, ul_map_y, rot_2, -gsd_y]
return ul_map_x, gsd_x, rot_1, ul_map_y, rot_2, -gsd_y
def to_mapinfo(self):
"""Return ENVI map info list using the attributes of the Geocoding instance.
......@@ -213,17 +220,75 @@ def geotransform2mapinfo(gt, prj):
:returns: ENVI geo info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
:rtype: list
"""
return Geocoding(gt=gt, prj=prj).to_mapinfo()
try:
return Geocoding(gt=gt, prj=prj).to_mapinfo()
except KeyError: # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
with NamedTemporaryFile(prefix='py_tools_ds__', suffix='.bsq') as tf_bsq:
path_hdr = os.path.splitext(tf_bsq.name)[0] + '.hdr'
try:
ds_inmem = get_GDAL_ds_inmem(np.array([[1, 2], [2, 3]]), gt, prj)
ds_out = gdal.GetDriverByName('ENVI').CreateCopy(tf_bsq.name, ds_inmem)
ds_out.FlushCache()
del ds_inmem, ds_out
with open(path_hdr, 'r') as inF:
content = inF.read()
if 'map info' in content:
res = re.search("map info = {(.*?)}", content, re.I).group(1)
map_info = [i.strip() for i in res.split(',')]
for i, ele in enumerate(map_info):
try:
map_info[i] = float(ele)
except ValueError:
pass
else:
map_info = ['Arbitrary', 1.0, 1.0, 0.0, 0.0, 1.0, 1.0]
finally:
os.remove(path_hdr)
return map_info
def mapinfo2geotransform(map_info):
# type: (Union[list, None]) -> list
# type: (Union[list, None]) -> tuple
"""Builds GDAL GeoTransform tuple from an ENVI geo info.
:param map_info: ENVI geo info (list), e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
:returns: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
"""
return Geocoding(mapinfo=map_info).to_geotransform()
try:
return Geocoding(mapinfo=map_info).to_geotransform()
except (KeyError, ValueError): # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
with NamedTemporaryFile(prefix='py_tools_ds__', suffix='.bsq') as tf_bsq:
path_hdr = os.path.splitext(tf_bsq.name)[0] + '.hdr'
try:
ds_out = gdal.GetDriverByName('ENVI').Create(tf_bsq.name, 2, 2, 1)
ds_out.GetRasterBand(1).WriteArray(np.array([[1, 2], [2, 3]]))
ds_out.FlushCache()
del ds_out
with open(path_hdr, 'r') as InHdr:
lines = InHdr.readlines()
lines.append('map info = { %s }\n' % ', '.join([str(i) for i in map_info]))
with open(path_hdr, 'w') as OutHdr:
OutHdr.writelines(lines)
ds = gdal.Open(tf_bsq.name)
gt = ds.GetGeoTransform()
del ds
return gt
finally:
os.remove(path_hdr)
def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
......
# -*- coding: utf-8 -*-
import numpy as np
import warnings
import multiprocessing
......@@ -20,9 +19,9 @@ from rasterio.warp import reproject as rio_reproject
from rasterio.warp import calculate_default_transform as rio_calc_transform
from rasterio.warp import Resampling
from pyresample.geometry import AreaDefinition, SwathDefinition
from pyresample.utils import get_area_def
from pyresample.bilinear import resample_bilinear
from pyresample.kd_tree import resample_nearest, resample_gauss, resample_custom
from pyresample.utils import get_area_def
from ...dtypes.conversion import dTypeDic_NumPy2GDAL
from ..projection import EPSG2WKT, WKT2EPSG, isProjectedOrGeographic, prj_equal, proj4_to_WKT
......@@ -545,6 +544,13 @@ class SensorMapGeometryTransformer(object):
del self.opts['radius_of_influence']
self.opts['radius'] = radius_of_influence
# NOTE: If pykdtree is built with OpenMP support (default) the number of threads is controlled with the
# standard OpenMP environment variable OMP_NUM_THREADS. The nprocs argument has no effect on pykdtree.
if 'nprocs' in self.opts:
if self.opts['nprocs'] > 1:
os.environ['OMP_NUM_THREADS'] = '%d' % opts['nprocs']
del self.opts['nprocs']
self.lats = lats
self.lons = lons
self.swath_definition = SwathDefinition(lons=lons, lats=lats)
......@@ -832,18 +838,10 @@ class SensorMapGeometryTransformer3D(object):
# define number of CPUs to use (but avoid sub-multiprocessing)
# -> parallelize either over bands or over image tiles
ncpus_avail = multiprocessing.cpu_count()
if 'nprocs' in opts and opts['nprocs']:
if self.lons.shape[2] > opts['nprocs']:
self.CPUs = opts['nprocs'] # parallelize over bands
opts['nprocs'] = 1
else:
self.CPUs = 1 # parallelize over imgage tiles
elif self.lons.shape[2] > ncpus_avail:
self.CPUs = ncpus_avail # parallelize over bands
else:
self.CPUs = 1
opts['nprocs'] = ncpus_avail # parallelize over imgage tiles
# bands: multiprocessing uses multiprocessing.Pool, implemented in to_map_geometry / to_sensor_geometry
# tiles: multiprocessing uses OpenMP implemented in pykdtree which is used by pyresample
self.opts['nprocs'] = opts.get('nprocs', multiprocessing.cpu_count())
self.mp_alg = 'bands' if self.lons.shape[2] >= opts['nprocs'] else 'tiles'
@staticmethod
def _to_map_geometry_2D(kwargs_dict):
......@@ -894,12 +892,16 @@ class SensorMapGeometryTransformer3D(object):
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_extent = tgt_extent or self._get_common_target_extent(tgt_epsg)
init_opts = self.opts.copy()
if self.mp_alg == 'bands':
del init_opts['nprocs'] # avoid sub-multiprocessing
args = [dict(
lons_2D=self.lons[:, :, band],
lats_2D=self.lats[:, :, band],
resamp_alg=self.resamp_alg,
radius_of_influence=self.radius_of_influence,
init_opts=self.opts,
init_opts=init_opts,
data_sensor_geo_2D=data[:, :, band],
tgt_prj=tgt_prj,
tgt_extent=tgt_extent,
......@@ -907,8 +909,8 @@ class SensorMapGeometryTransformer3D(object):
band_idx=band
) for band in range(data.shape[2])]
if self.CPUs > 1:
with multiprocessing.Pool(self.CPUs) as pool:
if self.mp_alg == 'bands':
with multiprocessing.Pool(self.opts['nprocs']) as pool:
result = pool.map(self._to_map_geometry_2D, args)
else:
result = [self._to_map_geometry_2D(argsdict) for argsdict in args]
......@@ -945,20 +947,24 @@ class SensorMapGeometryTransformer3D(object):
if data.ndim != 3:
raise ValueError(data.ndim, "'data' must have 3 dimensions.")
init_opts = self.opts.copy()
if self.mp_alg == 'bands':
del init_opts['nprocs'] # avoid sub-multiprocessing
args = [dict(
lons_2D=self.lons[:, :, band],
lats_2D=self.lats[:, :, band],
resamp_alg=self.resamp_alg,
radius_of_influence=self.radius_of_influence,
init_opts=self.opts,
init_opts=init_opts,
data_map_geo_2D=data[:, :, band],
src_prj=src_prj,
src_extent=src_extent,
band_idx=band
) for band in range(data.shape[2])]
if self.CPUs > 1:
with multiprocessing.Pool(self.CPUs) as pool:
if self.mp_alg == 'bands':
with multiprocessing.Pool(self.opts['nprocs']) as pool:
result = pool.map(self._to_sensor_geometry_2D, args)
else:
result = [self._to_sensor_geometry_2D(argsdict) for argsdict in args]
......
__version__ = '0.14.9'
__versionalias__ = '20190214_01'
__version__ = '0.14.10'
__versionalias__ = '20190219_01'
......@@ -11,13 +11,14 @@ Tests for `py_tools_ds.geo.map_info` module.
import unittest
from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
from py_tools_ds.geo.projection import EPSG2WKT
geotransform_utm = [331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0]
geotransform_utm_rotated = [331185.0, 12.202099292274006, 27.406363729278027,
5840115.0, 27.406363729278027, -12.202099292274006]
geotransform_local = [0, 1, 0, 0, 0, -1]
geotransform_local_rotated = [0.0, 6.123233995736766e-17, 1.0, 0.0, 1.0, -6.123233995736766e-17]
geotransform_utm = (331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0)
geotransform_utm_rotated = (331185.0, 12.202099292274006, 27.406363729278027,
5840115.0, 27.406363729278027, -12.202099292274006)
geotransform_local = (0, 1, 0, 0, 0, -1)
geotransform_local_rotated = (0.0, 6.123233995736766e-17, 1.0, 0.0, 1.0, -6.123233995736766e-17)
map_info_utm = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84']
map_info_utm_rotated = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84', 'rotation=66.00000']
map_info_local = ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North']
......@@ -77,6 +78,18 @@ class Test_geotransform2mapinfo(unittest.TestCase):
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, exp_map_info)
def test_prj_is_not_wgs84(self):
exp_map_info = ['Lambert Azimuthal Equal Area', 1.0, 1.0, 4526026.0, 3284919.5, 10.0, 10.0]
map_info = geotransform2mapinfo(gt=(4526026.0, 10.0, 0.0, 3284919.5, 0.0, -10.0), prj=EPSG2WKT(3035))
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, exp_map_info)
def test_prj_is_not_wgs84_rotation(self):
exp_map_info = ['Lambert Azimuthal Equal Area', 1.0, 1.0, 331185.0, 5840115.0, 30.0, 30.0, 'rotation=66']
map_info = geotransform2mapinfo(gt=geotransform_utm_rotated, prj=EPSG2WKT(3035))
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, exp_map_info)
def test_gt_contains_rotation(self):
map_info = geotransform2mapinfo(gt=geotransform_utm_rotated, prj=wkt_utm)
self.assertTrue(isinstance(map_info, list))
......@@ -92,25 +105,30 @@ class Test_mapinfo2geotransform(unittest.TestCase):
def test_map_info_is_valid(self):
gt = mapinfo2geotransform(map_info_utm)
self.assertTrue(isinstance(gt, list))
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, geotransform_utm)
# test
gt = mapinfo2geotransform(['Arbitrary', 1, 1, 5, -7, 1, 1, 0, 'North'])
self.assertTrue(isinstance(gt, list))
self.assertEqual(gt, [5, 1, 0, -7, 0, -1])
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, (5, 1, 0, -7, 0, -1))
def test_map_info_is_empty(self):
gt = mapinfo2geotransform(None)
self.assertTrue(isinstance(gt, list))
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, geotransform_local)
def test_map_info_contains_LAEA_proj(self):
gt = mapinfo2geotransform(['Lambert Azimuthal Equal Area', 1.0, 1.0, 4526026.0, 3284919.5, 10.0, 10.0])
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, (4526026.0, 10.0, 0.0, 3284919.5, 0.0, -10.0))
def test_map_info_contains_rotation(self):
gt = mapinfo2geotransform(map_info_utm_rotated)
self.assertTrue(isinstance(gt, list))
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, geotransform_utm_rotated)
def test_map_info_is_local_contains_rotation(self):
gt = mapinfo2geotransform(map_info_local_rotated)
self.assertTrue(isinstance(gt, list))
self.assertTrue(isinstance(gt, (tuple, list)))
self.assertEqual(gt, geotransform_local_rotated)
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