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

Completed implementation of geo.map_info.Geocoding class. Added tests for rotated datasets.

parent ee68f68a
# -*- coding: utf-8 -*-
import re
import math
import warnings
from typing import Union # noqa F401 # flake8 issue
try:
......@@ -9,82 +10,21 @@ except ImportError:
# noinspection PyPackageRequirements
import osr
from .coord_trafo import transform_utm_to_wgs84, transform_any_prj
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).
: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
"""
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.')
UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
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
def utm2wgs84(utmX, utmY): return transform_utm_to_wgs84(utmX, utmY, srs.GetUTMZone())
def is_UTM_North_South(LonLat): return 'North' if LonLat[1] >= 0. else 'South'
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]
else:
map_info = ['Arbitrary', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), 0, 'North']
return map_info
def mapinfo2geotransform(map_info):
# type: (Union[list, None]) -> list
"""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]
"""
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]
class Geocoding(object):
def __init__(self):
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
"""
self.prj_name = 'Arbitrary'
self.ul_x_map = 0.
self.ul_y_map = 0.
......@@ -101,8 +41,34 @@ class Geocoding(object):
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])
: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)):
......@@ -111,15 +77,26 @@ class Geocoding(object):
raise ValueError("'gt' must contain 6 elements.")
self.ul_x_map = float(gt[0])
self.gsd_x = float(gt[1])
self.ul_y_map = float(gt[3])
self.gsd_y = float(abs(gt[5]))
# handle rotations
self.rot_1_rad = math.asin(gt[2] / self.gsd_x)
self.rot_1_deg = math.degrees(self.rot_1_deg)
self.rot_2_rad = math.asin(gt[4] / self.gsd_y)
self.rot_2_deg = math.degrees(self.rot_2_rad)
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)
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()
......@@ -146,6 +123,12 @@ class Geocoding(object):
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
......@@ -182,6 +165,11 @@ class Geocoding(object):
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))
......@@ -195,6 +183,10 @@ class Geocoding(object):
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
......@@ -213,11 +205,24 @@ class Geocoding(object):
def geotransform2mapinfo(gt, prj):
return Geocoding().from_geotransform_projection(gt, prj).to_mapinfo()
# 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):
return Geocoding().from_mapinfo(map_info).to_geotransform()
# type: (Union[list, None]) -> list
"""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()
def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
......
......@@ -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 = \
"""
......@@ -72,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):
......@@ -89,3 +104,14 @@ 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)
print(gt)
self.assertTrue(isinstance(gt, list))
self.assertEqual(gt, geotransform_local_rotated)
Supports Markdown
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