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

Daniel Scheffler's avatar
Daniel Scheffler committed
3
from typing import Union  # noqa F401  # flake8 issue
Daniel Scheffler's avatar
Daniel Scheffler committed
4
5
6
try:
    from osgeo import osr
except ImportError:
Daniel Scheffler's avatar
Daniel Scheffler committed
7
    # noinspection PyPackageRequirements
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
    import osr

from .coord_trafo import transform_utm_to_wgs84

12
13
14
15
__author__ = "Daniel Scheffler"


def geotransform2mapinfo(gt, prj):
16
    """Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
Daniel Scheffler's avatar
Daniel Scheffler committed
17
18
    :param gt:  GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
    :param prj: GDAL Projection - WKT Format
19
    :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
20
21
22
    :rtype:     list
    """

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

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

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

    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
45
46
47
48
    return map_info


def mapinfo2geotransform(map_info):
Daniel Scheffler's avatar
Daniel Scheffler committed
49
    # type: (list) -> Union[list, None]
50
    """Builds GDAL GeoTransform tuple from an ENVI geo info.
51
52

    :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
53
54
    :returns:        GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
    """
55
56
    if map_info:
        # FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
57
58
59
60
61
        # 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)

62
63
64
65
66
        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]))
67

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

Daniel Scheffler's avatar
Daniel Scheffler committed
70
71
72
73
74

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'."
75

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

Daniel Scheffler's avatar
Daniel Scheffler committed
81
82
    for px in xarr:
        for py in yarr:
83
84
85
            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
86
        yarr.reverse()
87
    del gdal_ds_GT
88

89
    return ext