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

Added support for rotated datasets: Implemented class...

Added support for rotated datasets: Implemented class geo.map_info.Geocoding(). Reimplemented geotransform2mapinfo() and mapinfo2geotransform(). Bugfix for geo.projection.isLocal()
parent b10b7fe2
# -*- coding: utf-8 -*-
import re
import math
from typing import Union # noqa F401 # flake8 issue
try:
......@@ -7,7 +9,8 @@ except ImportError:
# noinspection PyPackageRequirements
import osr
from .coord_trafo import transform_utm_to_wgs84
from .coord_trafo import transform_utm_to_wgs84, transform_any_prj
from .projection import get_proj4info, proj4_to_dict, isLocal
__author__ = "Daniel Scheffler"
......@@ -80,6 +83,143 @@ def mapinfo2geotransform(map_info):
return [0, 1, 0, 0, 0, -1]
class Geocoding(object):
def __init__(self):
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 = ''
def from_geotransform_projection(self, gt, prj):
# type: (Union[list, tuple], str) -> self
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.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)
# handle projection
srs = osr.SpatialReference()
srs.ImportFromWkt(prj)
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
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'
del srs
return self
def from_mapinfo(self, mapinfo):
# 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:
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):
# 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):
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])
# 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):
return Geocoding().from_geotransform_projection(gt, prj).to_mapinfo()
def mapinfo2geotransform(map_info):
return Geocoding().from_mapinfo(map_info).to_geotransform()
def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
"""Returns (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform."""
assert gdal_ds or (gt and cols and rows), \
......
......@@ -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()
......
......@@ -54,12 +54,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))
......
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