map_info.py 13.4 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
from tempfile import TemporaryDirectory
30

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

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

44
45
46
__author__ = "Daniel Scheffler"


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

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

80
    def from_geotransform_projection(self, gt, prj):
81
        # type: (Union[list, tuple], str) -> 'Geocoding'
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
        """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
        """
101
102
103
104
105
106
107
108
109
110
111
        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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
            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)
127
128
129
130
131
132
133
134

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

            if isLocal(prj):
                self.prj_name = 'Arbitrary'
            else:
135
136
137
138
139
140
141
142
143
                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?
144
                    #       -> use CRS.datum.name instead?
145
146
                    self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum']  # proj4['ellps']?
                    self.units = proj4['unit'] if 'unit' in proj4 else self.units
147
148
149
150
151
152
153
154
155
156
157

                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):
158
        # type: (Union[list, tuple]) -> 'Geocoding'
159
160
161
162
163
        """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
        """
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
197
198
        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):
199
        # type: () -> tuple
200
201
        """Return GDAL GeoTransform list using the attributes of the Geocoding instance.

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

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

217
        return ul_map_x, gsd_x, rot_1, ul_map_y, rot_2, -gsd_y
218
219

    def to_mapinfo(self):
220
221
222
223
        """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 ]
        """
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
        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):
242
243
244
245
246
247
248
    # 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
    """
249
250
251
252
    try:
        return Geocoding(gt=gt, prj=prj).to_mapinfo()

    except KeyError:  # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
253
254
255
        if int(gdal.__version__[0]) < 3:
            # noinspection PyTypeChecker
            prj = CRS(prj).to_wkt(version="WKT1_GDAL")
256

257
258
259
260
        with TemporaryDirectory() as fdir:
            fn = "py_tools_ds__geotransform2mapinfo_temp"
            fP_bsq = os.path.join(fdir, fn + '.bsq')
            fP_hdr = os.path.join(fdir, fn + '.hdr')
261

262
            ds_out = gdal.GetDriverByName('ENVI').Create(fP_bsq, 2, 2, 1, gdal.GDT_Int32)
263
264
265
            ds_out.SetGeoTransform(gt)
            ds_out.SetProjection(prj)
            ds_out.FlushCache()
266
            del ds_out
267

268
            with open(fP_hdr, 'r') as inF:
269
270
271
272
273
274
275
276
277
278
279
280
281
282
                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]

        return map_info
283
284
285


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

    except (KeyError, ValueError):  # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection
296
297
298
299
        with TemporaryDirectory() as fdir:
            fn = "py_tools_ds__mapinfo2geotransform_temp"
            fP_bsq = os.path.join(fdir, fn + '.bsq')
            fP_hdr = os.path.join(fdir, fn + '.hdr')
300

301
            ds_out = gdal.GetDriverByName('ENVI').Create(fP_bsq, 2, 2, 1)
302
303
            ds_out.GetRasterBand(1).WriteArray(np.array([[1, 2], [2, 3]]))
            ds_out.FlushCache()
304
            del ds_out
305

306
            with open(fP_hdr, 'r') as InHdr:
307
                lines = InHdr.readlines()
308

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

311
            with open(fP_hdr, 'w') as OutHdr:
312
                OutHdr.writelines(lines)
313

314
            ds = gdal.Open(fP_bsq)
315
316
            gt = ds.GetGeoTransform()
            del ds
317

318
        return gt
319
320


Daniel Scheffler's avatar
Daniel Scheffler committed
321
322
323
324
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'."
325

Daniel Scheffler's avatar
Daniel Scheffler committed
326
    gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
327
328
329
    ext = []
    xarr = [0, gdal_ds.RasterXSize if gdal_ds else cols]
    yarr = [0, gdal_ds.RasterYSize if gdal_ds else rows]
330

Daniel Scheffler's avatar
Daniel Scheffler committed
331
332
    for px in xarr:
        for py in yarr:
333
334
335
            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
336
        yarr.reverse()
337
    del gdal_ds_GT
338

339
    return ext