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

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/>.

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

30
31
__author__ = "Daniel Scheffler"

Daniel Scheffler's avatar
Daniel Scheffler committed
32
33

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


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


Daniel Scheffler's avatar
Daniel Scheffler committed
49
def transform_wgs84_to_utm(lon, lat, zone=None):
50
51
52
53
54
    # type: (float, float, int) -> (float, float)
    """Return easting, northing, altitude for the given Lon/Lat coordinate.
    :param lon:     longitude
    :param lat:     latitude
    :param zone:    UTM zone (if not set it is derived automatically)
Daniel Scheffler's avatar
Daniel Scheffler committed
55
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
56
57
58
    if not (-180. <= lon <= 180. and -90. <= lat <= 90.):
        raise ValueError((lon, lat), 'Input coordinates for transform_wgs84_to_utm() out of range.')

59
60
61
    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
62
63
64
    return UTM(lon, lat)


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

    :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>)
71
72
    :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>)
73
74
    :return:
    """
75
    transformer = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(prj_src),
76
77
                                              pyproj.CRS.from_user_input(prj_tgt),
                                              always_xy=True)
78
    return transformer.transform(x, y)
Daniel Scheffler's avatar
Daniel Scheffler committed
79
80


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

85
86
87
88
89
90
91
92
93
    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:
    """
94
95
    drv = gdal.GetDriverByName('MEM')
    geoloc_ds = drv.Create('geoloc', Xarr.shape[1], Xarr.shape[0], 3, gdal.GDT_Float64)
96
97
98
99
100
    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
101
102
103
104
105
106
107
108
109
110
    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),
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
                                               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
126
def mapXY2imXY(mapXY, gt):
127
128
    # 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
129

130
    :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
131
    :param gt:      <list> GDAL geotransform
132
    :returns:       <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2]
Daniel Scheffler's avatar
Daniel Scheffler committed
133
    """
134
    if isinstance(mapXY, np.ndarray):
135
136
        ndim = mapXY.ndim
        if ndim == 1:
137
138
            mapXY = mapXY.reshape(1, 2)
        assert mapXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % mapXY.shape
139
        imXY = np.empty_like(mapXY, dtype=np.float)
140
141
        imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1]
        imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5]
142
        return imXY if ndim > 1 else imXY.flatten()
143
    else:
144
145
        return (mapXY[0] - gt[0]) / gt[1],\
               (mapXY[1] - gt[3]) / gt[5]
Daniel Scheffler's avatar
Daniel Scheffler committed
146
147
148


def imXY2mapXY(imXY, gt):
149
150
    # 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
151

152
    :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
153
    :param gt:      <list> GDAL geotransform
154
    :returns:       <tuple, np.ndarray> geo coordinate tuple X/Y (mapX, mapY) or np.ndarray [Nx2]
Daniel Scheffler's avatar
Daniel Scheffler committed
155
    """
156
    if isinstance(imXY, np.ndarray):
157
        ndim = imXY.ndim
158
159
160
        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
161
162
163
164
        mapXY = np.empty_like(imXY, dtype=np.float)
        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()
165
    else:
166
167
        return (gt[0] + imXY[0] * abs(gt[1])),\
               (gt[3] - imXY[1] * abs(gt[5]))
Daniel Scheffler's avatar
Daniel Scheffler committed
168
169


170
def mapYX2imYX(mapYX, gt):
171
172
    return (mapYX[0] - gt[3]) / gt[5],\
           (mapYX[1] - gt[0]) / gt[1]
173
174
175


def imYX2mapYX(imYX, gt):
176
177
    return gt[3] - (imYX[0] * abs(gt[5])),\
           gt[0] + (imYX[1] * gt[1])
178
179


180
181
182
183
def pixelToMapYX(pixelCoords, geotransform, projection):
    # MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform)

    # Create a spatial reference object
184
    srs = osr.SpatialReference()
185
186
    srs.ImportFromWkt(projection)

187
    # Set up the coordinate transformation object
188
    ct = osr.CoordinateTransformation(srs, srs)
189

190
    pixelCoords = [pixelCoords] if not type(pixelCoords[0]) in [list, tuple] else pixelCoords
191

192
193
194
    mapXYPairs = [ct.TransformPoint(pixX * geotransform[1] + geotransform[0],
                                    pixY * geotransform[5] + geotransform[3])[:2]
                  for pixX, pixY in pixelCoords]
195

196
    return [[mapY, mapX] for mapX, mapY in mapXYPairs]
197
198


199
200
201
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
202
203
204

    :param pixelPairs:      Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]]
    :param geotransform:    GDAL GeoTransform
205
    :param projection:      WKT string
Daniel Scheffler's avatar
Daniel Scheffler committed
206
207
    :returns:               The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]]
    """
208
209
210
211
212
213
214
215
216
217
218
219
    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
220

221
222
223
224

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.
225
226

    :param latLonPairs:     The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]]
Daniel Scheffler's avatar
Daniel Scheffler committed
227
    :param geotransform:    GDAL GeoTransform
228
    :param projection:      WKT string
Daniel Scheffler's avatar
Daniel Scheffler committed
229
230
    :return:            The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]]
    """
231
232
233
234
235
236
237
238
239
240
241
242
243
    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
244
245
    return pixelPairs

246

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

    # apply inverse geo transform
253
254
255
    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
256

257
258
259
    return pixel_x, abs(pixel_y)  # y pixel is likly a negative value given geo_transform


Daniel Scheffler's avatar
Daniel Scheffler committed
260
261
262
263
264
265
266
267
268
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:
    """
269
270
    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
271
272
273
            for GCP in gcpList]


274
def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
275
    """Reproject any shapely geometry from one projection to another.
276
277
278
279
280

    :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>)
    """
281
282
    project = pyproj.Transformer.from_proj(
        pyproj.CRS.from_user_input(prj_src),
283
284
        pyproj.CRS.from_user_input(prj_tgt),
        always_xy=True)
285

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