Commit 7c187233 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added workaround for ETRS/LAEA projection incompatibility + tests.

parent 2021ae17
Pipeline #3721 failed with stages
in 1 minute and 25 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.
......@@ -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
......@@ -162,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:
......@@ -180,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.
......@@ -212,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):
......
......@@ -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