Skip to content
Snippets Groups Projects
Commit 7447bfa8 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'feature/add_rotation_support'

parents b10b7fe2 168e3dd1
Branches
Tags
No related merge requests found
Pipeline #
Show whitespace changes
Inline Side-by-side
py_tools_ds/geo/map_info.py
+ 195
− 50
View file @ 7447bfa8
# -*- coding: utf-8 -*-
import re
import math
import warnings
from typing import Union # noqa F401 # flake8 issue
try:
......@@ -7,50 +10,209 @@ except ImportError:
# noinspection PyPackageRequirements
import osr
from .coord_trafo import transform_utm_to_wgs84
from .coord_trafo import transform_any_prj
from .projection import get_proj4info, proj4_to_dict, isLocal
__author__ = "Daniel Scheffler"
def geotransform2mapinfo(gt, prj):
# type: (Union[list, tuple, None], Union[str, None]) -> list
"""Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
class Geocoding(object):
def __init__(self, mapinfo=None, gt=None, prj=''):
# type: (Union[list, tuple], Union[list, tuple], str) -> None
"""Get an instance of the Geocoding object.
:param mapinfo: ENVI map info, e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
:param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
:param prj: GDAL Projection - WKT Format
:returns: ENVI geo info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
:rtype: list
"""
self.prj_name = 'Arbitrary'
self.ul_x_map = 0.
self.ul_y_map = 0.
self.ul_x_px = 1.
self.ul_y_px = 1.
self.gsd_x = 1.
self.gsd_y = 1.
self.rot_1_deg = 0.
self.rot_1_rad = 0.
self.rot_2_deg = 0.
self.rot_2_rad = 0.
self.utm_zone = 0
self.utm_north_south = 'North'
self.datum = ''
self.units = ''
if mapinfo:
if gt or prj:
warnings.warn("'gt' and 'prj' are not respected if mapinfo is given!")
self.from_mapinfo(mapinfo)
elif gt:
self.from_geotransform_projection(gt, prj)
def from_geotransform_projection(self, gt, prj):
# type: (Union[list, tuple], str) -> self
"""Create Geocoding object from GDAL GeoTransform + WKT projection string.
HOW COMPUTATION OF RADIANTS WORKS:
Building on top of the computation within self.to_geotransform():
gt[1] = math.cos(rotation_rad) * gsd_x
gt[2] = math.sin(rotation_rad) * gsd_x
-> we have to solve this equation system to get rotation_rad:
gsd_x = gt[2] / math.sin(rotation_rad)
gt[1] = math.cos(rotation_rad) * gt[2] / math.sin(rotation_rad)
gt[1] * math.sin(rotation_rad) = math.cos(rotation_rad) * gt[2]
math.sin(rotation_rad) / math.cos(rotation_rad) = gt[2] / gt[1]
math.tan(rotation_rad) = gt[2] / gt[1]
rotation_rad = math.atan(gt[2] / gt[1])
if gt in [None, [0, 1, 0, 0, 0, -1], (0, 1, 0, 0, 0, -1)]:
return ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North']
if gt[2] != 0 or gt[4] != 0: # TODO
raise NotImplementedError('Currently rotated datasets are not supported.')
:param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
:param prj: GDAL Projection - WKT Format
:return: instance of Geocoding
"""
if gt not in [None, [0, 1, 0, 0, 0, -1], (0, 1, 0, 0, 0, -1)]:
# validate input geotransform
if not isinstance(gt, (list, tuple)):
raise TypeError("'gt' must be a list or a tuple. Received type %s." % type(gt))
if len(gt) != 6:
raise ValueError("'gt' must contain 6 elements.")
self.ul_x_map = float(gt[0])
self.ul_y_map = float(gt[3])
self.gsd_y = float(abs(gt[5]))
# handle rotations
if float(gt[2]) == 0:
# no rotation. use default angles from init
self.gsd_x = float(gt[1])
else:
self.rot_1_rad = math.atan(gt[2] / gt[1])
self.rot_1_deg = math.degrees(self.rot_1_rad)
self.gsd_x = gt[2] / math.sin(self.rot_1_rad)
UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
if float(gt[4]) == 0:
# no rotation. use default angles from init
self.gsd_y = float(abs(gt[5]))
else:
self.rot_2_rad = math.atan(gt[4] / gt[5])
self.rot_2_deg = math.degrees(self.rot_2_rad)
self.gsd_y = gt[4] / math.sin(self.rot_2_rad)
# handle projection
srs = osr.SpatialReference()
srs.ImportFromWkt(prj)
if prj and not srs.IsLocal():
Proj4 = [i[1:] for i in srs.ExportToProj4().split()]
Proj4_proj = [v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] == 'proj'][0]
Proj4_ellps = \
[v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] in ['ellps', 'datum']][0]
proj = 'Geographic Lat/Lon' if Proj4_proj == 'longlat' else 'UTM' if Proj4_proj == 'utm' else Proj4_proj
ellps = 'WGS-84' if Proj4_ellps == 'WGS84' else Proj4_ellps
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']
self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum'] # proj4['ellps']?
self.units = proj4['unit'] if 'unit' in proj4 else self.units
def utm2wgs84(utmX, utmY): return transform_utm_to_wgs84(utmX, utmY, srs.GetUTMZone())
if self.prj_name == 'UTM':
self.utm_zone = srs.GetUTMZone()
self.utm_north_south = \
'North' if transform_any_prj(prj, 4326, self.ul_x_map, self.ul_y_map)[0] >= 0. else 'South'
def is_UTM_North_South(LonLat): return 'North' if LonLat[1] >= 0. else 'South'
del srs
north_south = is_UTM_North_South(utm2wgs84(UL_X, UL_Y))
map_info = [proj, 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), ellps] if proj != 'UTM' else \
['UTM', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), srs.GetUTMZone(), north_south, ellps]
return self
def from_mapinfo(self, mapinfo):
# type: (Union[list, tuple]) -> self
"""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)):
raise TypeError("'mapinfo' must be a list or a tuple. Received type %s." % type(mapinfo))
def assert_mapinfo_length(min_len):
if len(mapinfo) < min_len:
raise ValueError("A map info of type '%s' must contain at least %s elements. Received %s."
% (mapinfo[0], min_len, len(mapinfo)))
assert_mapinfo_length(10 if mapinfo[0] == 'UTM' else 9 if mapinfo[0] == 'Arbitrary' else 8)
# parse mapinfo
self.prj_name = mapinfo[0]
self.ul_x_px, self.ul_y_px, self.ul_x_map, self.ul_y_map, self.gsd_x = (float(i) for i in mapinfo[1:6])
self.gsd_y = float(abs(mapinfo[6]))
if self.prj_name == 'UTM':
self.utm_zone = mapinfo[7]
self.utm_north_south = mapinfo[8]
self.datum = mapinfo[9]
else:
map_info = ['Arbitrary', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), 0, 'North']
self.datum = mapinfo[7]
# handle rotation
for i in mapinfo:
if isinstance(i, str) and re.search('rotation', i, re.I):
self.rot_1_deg = float(i.split('=')[1].strip())
self.rot_2_deg = self.rot_1_deg
self.rot_1_rad = math.radians(self.rot_1_deg)
self.rot_2_rad = math.radians(self.rot_2_deg)
return self
def to_geotransform(self):
# type: () -> list
"""Return GDAL GeoTransform list using the attributes of the Geocoding instance.
:return: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
"""
# handle pixel coordinates of UL unequal to (1/1)
ul_map_x = self.ul_x_map if self.ul_x_px == 1 else (self.ul_x_map - (self.ul_x_px * self.gsd_x - self.gsd_x))
ul_map_y = self.ul_y_map if self.ul_y_px == 1 else (self.ul_y_map - (self.ul_y_px * self.gsd_y - self.gsd_y))
# handle rotation and pixel sizes
gsd_x, rot_1 = (self.gsd_x, 0) if self.rot_1_deg == 0 else \
(math.cos(self.rot_1_rad) * self.gsd_x, math.sin(self.rot_1_rad) * self.gsd_x)
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]
def to_mapinfo(self):
"""Return ENVI map info list using the attributes of the Geocoding instance.
:return: ENVI map info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
"""
mapinfo = [self.prj_name, self.ul_x_px, self.ul_y_px, self.ul_x_map, self.ul_y_map, self.gsd_x, abs(self.gsd_y)]
# add UTM infos
if self.prj_name in ['UTM', 'Arbitrary']:
mapinfo.extend([self.utm_zone, self.utm_north_south])
return map_info
# add datum
if self.prj_name != 'Arbitrary':
mapinfo.append(self.datum)
# add rotation
if self.rot_1_deg != 0.:
mapinfo.append('rotation=%.5f' % self.rot_1_deg)
return mapinfo
def geotransform2mapinfo(gt, prj):
# type: (Union[list, tuple, None], Union[str, None]) -> list
"""Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
:param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
:param prj: GDAL Projection - WKT Format
: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()
def mapinfo2geotransform(map_info):
......@@ -60,24 +222,7 @@ def mapinfo2geotransform(map_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]
"""
if map_info:
# FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
# FIXME in case of rotated datasets: map info contains additional key: {... , rotation=1.0} --> ANGLE IN DEGREE
# validate input map info
exp_len = 10 if map_info[0] == 'UTM' else 9 if map_info[0] == 'Arbitrary' else 8
assert isinstance(map_info, list) and len(map_info) == exp_len, \
"The map info argument has to be a list of length %s. Got %s." % (exp_len, len(map_info))
if map_info[1] == 1 and map_info[2] == 1:
ULmapX, ULmapY = float(map_info[3]), float(map_info[4])
else:
ULmapX = float(map_info[3]) - (float(map_info[1]) * float(map_info[5]) - float(map_info[5]))
ULmapY = float(map_info[4]) + (float(map_info[2]) * float(map_info[6]) - float(map_info[6]))
return [ULmapX, float(map_info[5]), 0., ULmapY, 0., -float(map_info[6])]
else:
return [0, 1, 0, 0, 0, -1]
return Geocoding(mapinfo=map_info).to_geotransform()
def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
......
py_tools_ds/geo/projection.py
+ 2
− 2
View file @ 7447bfa8
......@@ -117,10 +117,10 @@ def isLocal(prj):
srs.ImportFromEPSG(int(prj.split(':')[1]))
elif prj.startswith('+proj='):
srs.ImportFromProj4(prj)
elif prj.startswith('GEOGCS') or prj.startswith('PROJCS') or prj.startswith('LOCAL_CS'):
elif 'GEOGCS' in prj or 'PROJCS' in prj or 'LOCAL_CS' in prj:
srs.ImportFromWkt(prj)
else:
raise RuntimeError('Unknown input projection.')
raise RuntimeError('Unknown input projection: \n%s' % prj)
return srs.IsLocal()
......
tests/test_geo/test_map_info.py
+ 28
− 2
View file @ 7447bfa8
......@@ -14,9 +14,14 @@ from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
geotransform_utm = [331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0]
map_info_utm = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84']
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']
map_info_local_rotated = ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North', 'rotation=90.00000']
wkt_utm = \
"""
......@@ -54,12 +59,13 @@ class Test_geotransform2mapinfo(unittest.TestCase):
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, map_info_utm)
def test_gt_is_arbitrary_or_none(self):
def test_gt_is_none(self):
# test gt=None
map_info = geotransform2mapinfo(gt=None, prj=wkt_utm)
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, map_info_local)
def test_gt_is_arbitrary(self):
# test gt=[0, 1, 0, 0, 0, -1]
map_info = geotransform2mapinfo(gt=geotransform_local, prj=wkt_utm)
self.assertTrue(isinstance(map_info, list))
......@@ -71,6 +77,16 @@ class Test_geotransform2mapinfo(unittest.TestCase):
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))
self.assertEqual(map_info, map_info_utm_rotated)
def test_gt_contains_rotation_prj_is_local(self):
map_info = geotransform2mapinfo(gt=geotransform_local_rotated, prj='')
self.assertTrue(isinstance(map_info, list))
self.assertEqual(map_info, map_info_local_rotated)
class Test_mapinfo2geotransform(unittest.TestCase):
......@@ -88,3 +104,13 @@ class Test_mapinfo2geotransform(unittest.TestCase):
gt = mapinfo2geotransform(None)
self.assertTrue(isinstance(gt, list))
self.assertEqual(gt, geotransform_local)
def test_map_info_contains_rotation(self):
gt = mapinfo2geotransform(map_info_utm_rotated)
self.assertTrue(isinstance(gt, 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.assertEqual(gt, geotransform_local_rotated)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment