map_info.py 3.79 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
# -*- coding: utf-8 -*-

3
from typing import TYPE_CHECKING, Union, TypeVar
Daniel Scheffler's avatar
Daniel Scheffler committed
4
5
6
7
8
9
10
try:
    from osgeo import osr
except ImportError:
    import osr

from .coord_trafo import transform_utm_to_wgs84

11
12
13
14
__author__ = "Daniel Scheffler"

if TYPE_CHECKING:
    _T_mapinfo_out = TypeVar(Union[list, None])
Daniel Scheffler's avatar
Daniel Scheffler committed
15

16
17

def geotransform2mapinfo(gt, prj):
18
    """Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
Daniel Scheffler's avatar
Daniel Scheffler committed
19
20
    :param gt:  GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
    :param prj: GDAL Projection - WKT Format
21
    :returns:   ENVI geo info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
Daniel Scheffler's avatar
Daniel Scheffler committed
22
23
24
    :rtype:     list
    """

Daniel Scheffler's avatar
Daniel Scheffler committed
25
    if gt in [None, [0, 1, 0, 0, 0, -1]]:
Daniel Scheffler's avatar
Daniel Scheffler committed
26
27
        return None

28
    if gt[2] != 0 or gt[4] != 0:  # TODO
Daniel Scheffler's avatar
Daniel Scheffler committed
29
        raise NotImplementedError('Currently rotated datasets are not supported.')
30

31
    srs = osr.SpatialReference()
Daniel Scheffler's avatar
Daniel Scheffler committed
32
    srs.ImportFromWkt(prj)
33
34
35
36
37
38
    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
Daniel Scheffler's avatar
Daniel Scheffler committed
39
    UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
40
41
42
43
44
45
46

    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'

    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(), is_UTM_North_South(utm2wgs84(UL_X, UL_Y)), ellps]
Daniel Scheffler's avatar
Daniel Scheffler committed
47
48
49
50
    return map_info


def mapinfo2geotransform(map_info):
51
    # type: (list) -> _T_mapinfo_out
52
    """Builds GDAL GeoTransform tuple from an ENVI geo info.
53
54

    :param map_info: ENVI geo info (list), e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
Daniel Scheffler's avatar
Daniel Scheffler committed
55
56
    :returns:        GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
    """
57
58
    if map_info:
        # FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
59
60
61
62
63
        # validate input map info
        exp_len = 10 if map_info[0] == 'UTM' 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." % (map_info, exp_len)

64
65
66
67
68
        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]))
69

70
        return [ULmapX, float(map_info[5]), 0., ULmapY, 0., -float(map_info[6])]
Daniel Scheffler's avatar
Daniel Scheffler committed
71

Daniel Scheffler's avatar
Daniel Scheffler committed
72
73
74
75
76

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), \
        "GEOP.get_corner_coordinates: Missing argument! Please provide either 'gdal_ds' or 'gt', 'cols' AND 'rows'."
77

Daniel Scheffler's avatar
Daniel Scheffler committed
78
    gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
79
80
81
    ext = []
    xarr = [0, gdal_ds.RasterXSize if gdal_ds else cols]
    yarr = [0, gdal_ds.RasterYSize if gdal_ds else rows]
82

Daniel Scheffler's avatar
Daniel Scheffler committed
83
84
    for px in xarr:
        for py in yarr:
85
86
87
            x = gdal_ds_GT[0] + (px * gdal_ds_GT[1]) + (py * gdal_ds_GT[2])
            y = gdal_ds_GT[3] + (px * gdal_ds_GT[4]) + (py * gdal_ds_GT[5])
            ext.append([x, y])
Daniel Scheffler's avatar
Daniel Scheffler committed
88
        yarr.reverse()
89
    del gdal_ds_GT
90

91
    return ext