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

3
4
# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
# operations when handling geospatial raster and vector data as well as projections.
5
#
6
7
8
9
# Copyright (C) 2016-2021
# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
#   Germany (https://www.gfz-potsdam.de/)
10
11
12
13
14
#
# 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).
#
15
16
17
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
26

Daniel Scheffler's avatar
Daniel Scheffler committed
27
import pyproj
28
import numpy as np
29
from shapely.ops import transform
30
from typing import Union, Sequence, List  # noqa F401  # flake8 issue
31
from osgeo import gdal, osr
Daniel Scheffler's avatar
Daniel Scheffler committed
32

33
34
__author__ = "Daniel Scheffler"

Daniel Scheffler's avatar
Daniel Scheffler committed
35
36

def transform_utm_to_wgs84(easting, northing, zone, south=False):
37
    # type: (float, float, int, bool) -> (float, float)
Daniel Scheffler's avatar
Daniel Scheffler committed
38
    """Transform an UTM coordinate to lon, lat, altitude."""
39
40
41
42
    UTM = pyproj.Proj(proj='utm',
                      zone=abs(zone),
                      ellps='WGS84',
                      south=(zone < 0 or south))
Daniel Scheffler's avatar
Daniel Scheffler committed
43
44
45
    return UTM(easting, northing, inverse=True)


Daniel Scheffler's avatar
Daniel Scheffler committed
46
def get_utm_zone(longitude):
47
    # type: (float) -> int
Daniel Scheffler's avatar
Daniel Scheffler committed
48
49
50
51
    """Get the UTM zone corresponding to the given longitude."""
    return int(1 + (longitude + 180.0) / 6.0)


Daniel Scheffler's avatar
Daniel Scheffler committed
52
def transform_wgs84_to_utm(lon, lat, zone=None):
53
54
    # type: (float, float, int) -> (float, float)
    """Return easting, northing, altitude for the given Lon/Lat coordinate.
Daniel Scheffler's avatar
Daniel Scheffler committed
55

56
57
58
    :param lon:     longitude
    :param lat:     latitude
    :param zone:    UTM zone (if not set it is derived automatically)
Daniel Scheffler's avatar
Daniel Scheffler committed
59
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
60
61
62
    if not (-180. <= lon <= 180. and -90. <= lat <= 90.):
        raise ValueError((lon, lat), 'Input coordinates for transform_wgs84_to_utm() out of range.')

63
64
65
    UTM = pyproj.Proj(proj='utm',
                      zone=get_utm_zone(lon) if zone is None else zone,
                      ellps='WGS84')
Daniel Scheffler's avatar
Daniel Scheffler committed
66
67
68
    return UTM(lon, lat)


69
def transform_any_prj(prj_src, prj_tgt, x, y):
70
    # type: (Union[str, int], Union[str, int], Union[float, np.ndarray], Union[float, np.ndarray]) -> (Union[float, np.ndarray], Union[float, np.ndarray])  # noqa
71
    """Transform X/Y data from any source projection to any target projection.
72
73
74

    :param prj_src:     GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
    :param prj_tgt:     GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
75
76
    :param x:           X-coordinate(s) to be transformed (either <float> or 1D-<np.ndarray>)
    :param y:           X-coordinate(s) to be transformed (either <float> or 1D-<np.ndarray>)
77
78
    :return:
    """
79
    transformer = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(prj_src),
80
81
                                              pyproj.CRS.from_user_input(prj_tgt),
                                              always_xy=True)
82
    return transformer.transform(x, y)
Daniel Scheffler's avatar
Daniel Scheffler committed
83
84


85
def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
86
87
88
    # type: (str, str, np.ndarray, np.ndarray, np.ndarray) -> Sequence[np.ndarray]
    """Transform a geolocation array from one projection into another.

89
90
91
92
93
94
95
96
97
    HINT: This function is faster than transform_any_prj but works only for geolocation arrays.

    :param prj_src:     WKT string
    :param prj_tgt:     WKT string
    :param Xarr:        np.ndarray of shape (rows,cols)
    :param Yarr:        np.ndarray of shape (rows,cols)
    :param Zarr:        np.ndarray of shape (rows,cols)
    :return:
    """
98
99
    drv = gdal.GetDriverByName('MEM')
    geoloc_ds = drv.Create('geoloc', Xarr.shape[1], Xarr.shape[0], 3, gdal.GDT_Float64)
100
101
102
103
104
    geoloc_ds.GetRasterBand(1).WriteArray(Xarr)
    geoloc_ds.GetRasterBand(2).WriteArray(Yarr)
    if Zarr is not None:
        geoloc_ds.GetRasterBand(3).WriteArray(Zarr)

Daniel Scheffler's avatar
Daniel Scheffler committed
105
106
107
108
109
110
111
112
113
114
    def strip_wkt(wkt):
        return wkt\
            .replace('\n', '')\
            .replace('\r', '')\
            .replace(' ', '')

    transformer = gdal.Transformer(None, None, ['SRC_SRS=' + strip_wkt(prj_src),
                                                'DST_SRS=' + strip_wkt(prj_tgt)])
    status = transformer.TransformGeolocations(geoloc_ds.GetRasterBand(1),
                                               geoloc_ds.GetRasterBand(2),
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
                                               geoloc_ds.GetRasterBand(3))

    if status:
        raise Exception('Error transforming coordinate array:  ' + gdal.GetLastErrorMsg())

    Xarr = geoloc_ds.GetRasterBand(1).ReadAsArray()
    Yarr = geoloc_ds.GetRasterBand(2).ReadAsArray()

    if Zarr is not None:
        Zarr = geoloc_ds.GetRasterBand(3).ReadAsArray()
        return Xarr, Yarr, Zarr
    else:
        return Xarr, Yarr


Daniel Scheffler's avatar
Daniel Scheffler committed
130
def mapXY2imXY(mapXY, gt):
131
132
    # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray]
    """Translate the given geo coordinates into pixel locations according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
133

134
    :param mapXY:   <tuple, np.ndarray> The geo coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
Daniel Scheffler's avatar
Daniel Scheffler committed
135
    :param gt:      <list> GDAL geotransform
136
    :returns:       <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2]
Daniel Scheffler's avatar
Daniel Scheffler committed
137
    """
138
    if isinstance(mapXY, np.ndarray):
139
140
        ndim = mapXY.ndim
        if ndim == 1:
141
142
            mapXY = mapXY.reshape(1, 2)
        assert mapXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % mapXY.shape
143
        imXY = np.empty_like(mapXY, dtype=float)
144
145
        imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1]
        imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5]
146
        return imXY if ndim > 1 else imXY.flatten()
147
    else:
148
149
        return (mapXY[0] - gt[0]) / gt[1],\
               (mapXY[1] - gt[3]) / gt[5]
Daniel Scheffler's avatar
Daniel Scheffler committed
150
151
152


def imXY2mapXY(imXY, gt):
153
154
    # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray]
    """Translate the given pixel X/Y locations into geo coordinates according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
155

156
    :param imXY:    <tuple, np.ndarray> The image coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
Daniel Scheffler's avatar
Daniel Scheffler committed
157
    :param gt:      <list> GDAL geotransform
158
    :returns:       <tuple, np.ndarray> geo coordinate tuple X/Y (mapX, mapY) or np.ndarray [Nx2]
Daniel Scheffler's avatar
Daniel Scheffler committed
159
    """
160
    if isinstance(imXY, np.ndarray):
161
        ndim = imXY.ndim
162
163
164
        if imXY.ndim == 1:
            imXY = imXY.reshape(1, 2)
        assert imXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % imXY.shape
165
        mapXY = np.empty_like(imXY, dtype=float)
166
167
168
        mapXY[:, 0] = gt[0] + imXY[:, 0] * abs(gt[1])
        mapXY[:, 1] = gt[3] - imXY[:, 1] * abs(gt[5])
        return mapXY if ndim > 1 else mapXY.flatten()
169
    else:
170
171
        return (gt[0] + imXY[0] * abs(gt[1])),\
               (gt[3] - imXY[1] * abs(gt[5]))
Daniel Scheffler's avatar
Daniel Scheffler committed
172
173


174
def mapYX2imYX(mapYX, gt):
175
176
    return (mapYX[0] - gt[3]) / gt[5],\
           (mapYX[1] - gt[0]) / gt[1]
177
178
179


def imYX2mapYX(imYX, gt):
180
181
    return gt[3] - (imYX[0] * abs(gt[5])),\
           gt[0] + (imYX[1] * gt[1])
182
183


184
185
186
187
def pixelToMapYX(pixelCoords, geotransform, projection):
    # MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform)

    # Create a spatial reference object
188
    srs = osr.SpatialReference()
189
190
    srs.ImportFromWkt(projection)

191
    # Set up the coordinate transformation object
192
    ct = osr.CoordinateTransformation(srs, srs)
193

194
    pixelCoords = [pixelCoords] if not type(pixelCoords[0]) in [list, tuple] else pixelCoords
195

196
197
198
    mapXYPairs = [ct.TransformPoint(pixX * geotransform[1] + geotransform[0],
                                    pixY * geotransform[5] + geotransform[3])[:2]
                  for pixX, pixY in pixelCoords]
199

200
    return [[mapY, mapX] for mapX, mapY in mapXYPairs]
201
202


203
204
205
def pixelToLatLon(pixelPairs, geotransform, projection):
    # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]]
    """Translate the given pixel X/Y locations into latitude/longitude locations.
Daniel Scheffler's avatar
Daniel Scheffler committed
206
207
208

    :param pixelPairs:      Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]]
    :param geotransform:    GDAL GeoTransform
209
    :param projection:      WKT string
Daniel Scheffler's avatar
Daniel Scheffler committed
210
211
    :returns:               The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]]
    """
212
213
214
215
216
217
218
219
220
221
222
223
    MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform)
    lons, lats = transform_any_prj(projection, 4326, MapXYPairs[:, 0], MapXYPairs[:, 1])

    # validate output lat/lon values
    if not (-180 <= np.min(lons) <= 180) or not (-180 <= np.max(lons) <= 180):
        raise RuntimeError('Output longitude values out of bounds.')
    if not (-90 <= np.min(lats) <= 90) or not (-90 <= np.max(lats) <= 90):
        raise RuntimeError('Output latitudes values out of bounds.')

    LatLonPairs = np.vstack([lats, lons]).T.tolist()

    return LatLonPairs
Daniel Scheffler's avatar
Daniel Scheffler committed
224

225
226
227
228

def latLonToPixel(latLonPairs, geotransform, projection):
    # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]]
    """Translate the given latitude/longitude pairs into pixel X/Y locations.
229
230

    :param latLonPairs:     The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]]
Daniel Scheffler's avatar
Daniel Scheffler committed
231
    :param geotransform:    GDAL GeoTransform
232
    :param projection:      WKT string
Daniel Scheffler's avatar
Daniel Scheffler committed
233
234
    :return:            The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]]
    """
235
236
237
238
239
240
241
242
243
244
245
246
247
    latLonArr = np.array(latLonPairs)

    # validate input lat/lon values
    if not (-180 <= np.min(latLonArr[:, 1]) <= 180) or not (-180 <= np.max(latLonArr[:, 1]) <= 180):
        raise ValueError('Longitude value out of bounds.')
    if not (-90 <= np.min(latLonArr[:, 0]) <= 90) or not (-90 <= np.max(latLonArr[:, 0]) <= 90):
        raise ValueError('Latitude value out of bounds.')

    mapXs, mapYs = transform_any_prj(4326, projection, latLonArr[:, 1], latLonArr[:, 0])
    imXYarr = mapXY2imXY(np.vstack([mapXs, mapYs]).T, geotransform)

    pixelPairs = imXYarr.tolist()

Daniel Scheffler's avatar
Daniel Scheffler committed
248
249
    return pixelPairs

250

251
def lonlat_to_pixel(lon, lat, inverse_geo_transform):
252
    """Translate the given lon, lat to the grid pixel coordinates in data array (zero start)."""
253
    # transform to utm
254
255
256
    utm_x, utm_y = transform_wgs84_to_utm(lon, lat)

    # apply inverse geo transform
257
258
259
    pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
    pixel_x = int(pixel_x) - 1  # adjust to 0 start for array
    pixel_y = int(pixel_y) - 1  # adjust to 0 start for array
260

261
262
263
    return pixel_x, abs(pixel_y)  # y pixel is likly a negative value given geo_transform


Daniel Scheffler's avatar
Daniel Scheffler committed
264
265
266
267
268
269
270
271
272
def transform_GCPlist(gcpList, prj_src, prj_tgt):
    """

    :param gcpList:     <list> list of ground control points in the output coordinate system
                                to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...]
    :param prj_src:     WKT string
    :param prj_tgt:     WKT string
    :return:
    """
273
274
    return [gdal.GCP(*(list(transform_any_prj(prj_src, prj_tgt, GCP.GCPX, GCP.GCPY)) +
                       [GCP.GCPZ, GCP.GCPPixel, GCP.GCPLine]))  # Python 2.7 compatible syntax
275
276
277
            for GCP in gcpList]


278
def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
279
    """Reproject any shapely geometry from one projection to another.
280
281
282
283
284

    :param shapelyGeometry: any shapely geometry instance
    :param prj_src:         GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
    :param prj_tgt:         GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
    """
285
286
    project = pyproj.Transformer.from_proj(
        pyproj.CRS.from_user_input(prj_src),
287
288
        pyproj.CRS.from_user_input(prj_tgt),
        always_xy=True)
289

290
    return transform(project.transform, shapelyGeometry)  # apply projection