map_info.py 13.6 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# py_tools_ds
#
# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

24
25
import re
import math
26
import warnings
27
import os
Daniel Scheffler's avatar
Daniel Scheffler committed
28
from typing import Union  # noqa F401  # flake8 issue
29

30
from pyproj import CRS
31
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
32
33
try:
    from osgeo import osr
34
    from osgeo import gdal
Daniel Scheffler's avatar
Daniel Scheffler committed
35
except ImportError:
Daniel Scheffler's avatar
Daniel Scheffler committed
36
    # noinspection PyPackageRequirements
Daniel Scheffler's avatar
Daniel Scheffler committed
37
    import osr
38
    import gdal
Daniel Scheffler's avatar
Daniel Scheffler committed
39

40
from .coord_trafo import transform_any_prj
41
from .projection import isLocal
Daniel Scheffler's avatar
Daniel Scheffler committed
42

43
44
45
__author__ = "Daniel Scheffler"


46
class Geocoding(object):
47
48
49
50
51
52
53
54
    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
        """
55
        # FIXME this class only supports UTM and Geographic Lat/Lon projections
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
        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 = ''

72
73
74
75
76
77
78
        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)

79
    def from_geotransform_projection(self, gt, prj):
80
        # type: (Union[list, tuple], str) -> 'Geocoding'
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
        """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
        """
100
101
102
103
104
105
106
107
108
109
110
        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])

            # handle rotations
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
            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)
126
127
128
129
130
131
132
133

            # handle projection
            srs = osr.SpatialReference()
            srs.ImportFromWkt(prj)

            if isLocal(prj):
                self.prj_name = 'Arbitrary'
            else:
134
135
136
137
138
139
140
141
142
143
144
                with warnings.catch_warnings():
                    warnings.filterwarnings("ignore", message="You will likely lose important projection information")

                    # get prj_name and datum
                    proj4 = CRS(prj).to_dict()  # FIXME avoid to convert to proj4
                    self.prj_name = \
                        'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \
                        'UTM' if proj4['proj'] == 'utm' else proj4['proj']
                    # FIXME there is no 'datum' key in case of, e.g., LAEA projection  # still true?
                    self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum']  # proj4['ellps']?
                    self.units = proj4['unit'] if 'unit' in proj4 else self.units
145
146
147
148
149
150
151
152
153
154
155

                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):
156
        # type: (Union[list, tuple]) -> 'Geocoding'
157
158
159
160
161
        """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
        """
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
        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):
197
        # type: () -> tuple
198
199
        """Return GDAL GeoTransform list using the attributes of the Geocoding instance.

Daniel Scheffler's avatar
Daniel Scheffler committed
200
201
202
        For equations, see:
         https://gis.stackexchange.com/questions/229952/rotate-envi-hyperspectral-imagery-with-gdal/229962

203
204
        :return:    GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
        """
205
206
207
208
209
210
211
212
213
214
        # 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)

215
        return ul_map_x, gsd_x, rot_1, ul_map_y, rot_2, -gsd_y
216
217

    def to_mapinfo(self):
218
219
220
221
        """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 ]
        """
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
        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):
240
241
242
243
244
245
246
    # 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
    """
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    try:
        return Geocoding(gt=gt, prj=prj).to_mapinfo()

    except KeyError:  # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
        fn_bsq = "py_tools_ds__geotransform2mapinfo_temp.bsq"
        fn_hdr = os.path.splitext(fn_bsq)[0] + '.hdr'
        fdir = os.path.join(os.path.abspath(os.curdir))

        try:
            ds_out = gdal.GetDriverByName('ENVI').Create(fn_bsq, 2, 2, 1, gdal.GDT_Int32)
            ds_out.SetGeoTransform(gt)
            ds_out.SetProjection(prj)
            ds_out.FlushCache()

            # noinspection PyUnusedLocal
            ds_out = None

            with open(fn_hdr, 'r') as inF:
                content = inF.read()
                if 'map info' in content:
                    res = re.search("map info = {(.*?)}", content, re.I).group(1)
                    map_info = [i.strip() for i in res.split(',')]

                    for i, ele in enumerate(map_info):
                        try:
                            map_info[i] = float(ele)
                        except ValueError:
                            pass
                else:
                    map_info = ['Arbitrary', 1.0, 1.0, 0.0, 0.0, 1.0, 1.0]

        finally:
            for fn in [fn_bsq, fn_hdr]:
                if os.path.exists(os.path.join(fdir, fn)):
                    gdal.Unlink(os.path.join(fdir, fn))

        return map_info
284
285
286


def mapinfo2geotransform(map_info):
287
    # type: (Union[list, None]) -> tuple
288
289
290
291
292
    """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]
    """
293
294
295
296
    try:
        return Geocoding(mapinfo=map_info).to_geotransform()

    except (KeyError, ValueError):  # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
297
298
299
        fn_bsq = "py_tools_ds__geotransform2mapinfo_temp.bsq"
        fn_hdr = os.path.splitext(fn_bsq)[0] + '.hdr'
        fdir = os.path.join(os.path.abspath(os.curdir))
300

301
302
303
304
        try:
            ds_out = gdal.GetDriverByName('ENVI').Create(fn_bsq, 2, 2, 1)
            ds_out.GetRasterBand(1).WriteArray(np.array([[1, 2], [2, 3]]))
            ds_out.FlushCache()
305
306
307

            # noinspection PyUnusedLocal
            ds_out = None
308

309
310
            with open(fn_hdr, 'r') as InHdr:
                lines = InHdr.readlines()
311

312
            lines.append('map info = { %s }\n' % ', '.join([str(i) for i in map_info]))
313

314
315
            with open(fn_hdr, 'w') as OutHdr:
                OutHdr.writelines(lines)
316

317
318
319
            ds = gdal.Open(fn_bsq)
            gt = ds.GetGeoTransform()
            del ds
320

321
            return gt
322

323
324
325
326
        finally:
            for fn in [fn_bsq, fn_hdr]:
                if os.path.exists(os.path.join(fdir, fn)):
                    gdal.Unlink(os.path.join(fdir, fn))
327
328


Daniel Scheffler's avatar
Daniel Scheffler committed
329
330
331
332
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'."
333

Daniel Scheffler's avatar
Daniel Scheffler committed
334
    gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
335
336
337
    ext = []
    xarr = [0, gdal_ds.RasterXSize if gdal_ds else cols]
    yarr = [0, gdal_ds.RasterYSize if gdal_ds else rows]
338

Daniel Scheffler's avatar
Daniel Scheffler committed
339
340
    for px in xarr:
        for py in yarr:
341
342
343
            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
344
        yarr.reverse()
345
    del gdal_ds_GT
346

347
    return ext