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

3
4
5
from functools import partial
import warnings

Daniel Scheffler's avatar
Daniel Scheffler committed
6
import pyproj
7
8
import numpy as np
from shapely.geometry import Polygon
9
10
from shapely.ops import transform
from typing import Union, TypeVar, TYPE_CHECKING
Daniel Scheffler's avatar
Daniel Scheffler committed
11
12
13
14
15
16
17
18

try:
    from osgeo import osr
    from osgeo import gdal
except ImportError:
    import osr
    import gdal

19
from .projection import get_proj4info, proj4_to_dict
Daniel Scheffler's avatar
Daniel Scheffler committed
20

21
22
23
24
__author__ = "Daniel Scheffler"

if TYPE_CHECKING:
    _T_xycoords = TypeVar(Union[tuple, np.ndarray])
Daniel Scheffler's avatar
Daniel Scheffler committed
25
26
27
28
29
30
31
32
33
34
35
36


def transform_utm_to_wgs84(easting, northing, zone, south=False):
    #    ''' returns lon, lat, altitude '''
    #    utm_coordinate_system = osr.SpatialReference()
    #    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
    #    is_northern = northingcalc_FullDataset_corner_positions > 0
    #    utm_coordinate_system.SetUTM(zone, is_northern)
    #
    #    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
    #
    #    # create transform component
37
    #    utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system)
Daniel Scheffler's avatar
Daniel Scheffler committed
38
39
    #    print(utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0))
    #    return utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0)  #[Lon,Lat,meters above sea level]
40
    UTM = pyproj.Proj(proj='utm', zone=abs(zone), ellps='WGS84', south=(zone < 0 or south))
Daniel Scheffler's avatar
Daniel Scheffler committed
41
42
43
44
45
46
47
48
49
    return UTM(easting, northing, inverse=True)


def transform_wgs84_to_utm(lon, lat, zone=None):
    """ returns easting, northing, altitude. If zone is not set it is derived automatically.
    :param lon:
    :param lat:
    :param zone:
    """
50
51
52
53
    assert -180. <= lon <= 180. and -90. <= lat <= 90., 'Input coordinates for transform_wgs84_to_utm() out of range. '\
                                                        'Got %s/%s.' % (lon, lat)

    def get_utm_zone(longitude): return int(1 + (longitude + 180.0) / 6.0)
Daniel Scheffler's avatar
Daniel Scheffler committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
    #    def is_northern(latitude):
    #        """
    #        Determines if given latitude is a northern for UTM
    #        """
    #        if (latitude < 0.0):
    #            return 0
    #        else:
    #            return 1
    #
    #    utm_coordinate_system = osr.SpatialReference()
    #    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
    #    utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))
    #
    #    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
    #
    #    # create transform component
70
    #    wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system)
Daniel Scheffler's avatar
Daniel Scheffler committed
71
72
73
74
75
    #    return wgs84_to_utm_geo_transform.TransformPoint(lon, lat, 0) # [RW,HW,meters above sea level]
    UTM = pyproj.Proj(proj='utm', zone=get_utm_zone(lon) if zone is None else zone, ellps='WGS84')
    return UTM(lon, lat)


76
def transform_any_prj(prj_src, prj_tgt, x, y):
77
78
79
80
81
82
83
84
    """Transforms X/Y data from any source projection to any target projection.

    :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>)
    :param x:
    :param y:
    :return:
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
85
86
    prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
    prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
87
88
    x, y = pyproj.transform(prj_src, prj_tgt, x, y)
    return x, y
Daniel Scheffler's avatar
Daniel Scheffler committed
89
90


91
92
93
94
95
96
97
98
99
100
101
def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
    """Transforms geolocation arrays from one projection into another.
    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:
    """
102
103
    drv = gdal.GetDriverByName('MEM')
    geoloc_ds = drv.Create('geoloc', Xarr.shape[1], Xarr.shape[0], 3, gdal.GDT_Float64)
104
105
106
107
108
109
    geoloc_ds.GetRasterBand(1).WriteArray(Xarr)
    geoloc_ds.GetRasterBand(2).WriteArray(Yarr)
    if Zarr is not None:
        geoloc_ds.GetRasterBand(3).WriteArray(Zarr)

    transformer = gdal.Transformer(None, None, ['SRC_SRS=' + prj_src, 'DST_SRS=' + prj_tgt])
110
    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
    # type: (tuple,list) -> _T_xycoords
128
    """Translates 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
137
        if mapXY.ndim == 1:
            mapXY = mapXY.reshape(1, 2)
        assert mapXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % mapXY.shape
138
        imXY = np.empty_like(mapXY)
139
140
        imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1]
        imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5]
141
142
        return imXY
    else:
143
        return (mapXY[0] - gt[0]) / gt[1], (mapXY[1] - gt[3]) / gt[5]
Daniel Scheffler's avatar
Daniel Scheffler committed
144
145
146


def imXY2mapXY(imXY, gt):
147
    # type: (tuple,list) -> _T_xycoords
148
    """Translates given pixel locations into geo coordinates according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
149

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


166
167
168
169
170
171
def mapYX2imYX(mapYX, gt):
    return (mapYX[0] - gt[3]) / gt[5], (mapYX[1] - gt[0]) / gt[1]


def imYX2mapYX(imYX, gt):
    return gt[3] - (imYX[0] * abs(gt[5])), gt[0] + (imYX[1] * gt[1])
172
173
174
175
176
177


def pixelToMapYX(pixelCoords, path_im=None, geotransform=None, projection=None):
    assert path_im or geotransform and projection, \
        "pixelToMapYX: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
    if geotransform and projection:
178
        gt, proj = geotransform, projection
179
    else:
180
        ds = gdal.Open(path_im)
181
182
183
184
        gt, proj = ds.GetGeoTransform(), ds.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(proj)
    # Set up the coordinate transformation object
185
    ct = osr.CoordinateTransformation(srs, srs)
186
187

    mapYmapXPairs = []
188
    pixelCoords = [pixelCoords] if not type(pixelCoords[0]) in [list, tuple] else pixelCoords
189
190
191
192
193
194
195
196
197
198
199
200
201

    for point in pixelCoords:
        # Translate the pixel pairs into untranslated points
        u_mapX = point[0] * gt[1] + gt[0]
        u_mapY = point[1] * gt[5] + gt[3]
        # Transform the points to the space
        (mapX, mapY, holder) = ct.TransformPoint(u_mapX, u_mapY)
        # Add the point to our return array
        mapYmapXPairs.append([mapY, mapX])

    return mapYmapXPairs


Daniel Scheffler's avatar
Daniel Scheffler committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def pixelToLatLon(pixelPairs, path_im=None, geotransform=None, projection=None):
    """The following method translates given pixel locations into latitude/longitude locations on a given GEOTIF
    This method does not take into account pixel size and assumes a high enough
    image resolution for pixel size to be insignificant

    :param pixelPairs:      Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]]
    :param path_im:         The file location of the input image
    :param geotransform:    GDAL GeoTransform
    :param projection:      GDAL Projection
    :returns:               The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]]
    """
    assert path_im or geotransform and projection, \
        "GEOP.pixelToLatLon: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
    if geotransform and projection:
216
        gt, proj = geotransform, projection
Daniel Scheffler's avatar
Daniel Scheffler committed
217
    else:
218
        ds = gdal.Open(path_im)
Daniel Scheffler's avatar
Daniel Scheffler committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
        gt, proj = ds.GetGeoTransform(), ds.GetProjection()

    # Create a spatial reference object for the dataset
    srs = osr.SpatialReference()
    srs.ImportFromWkt(proj)
    # Set up the coordinate transformation object
    srsLatLong = srs.CloneGeogCS()
    ct = osr.CoordinateTransformation(srs, srsLatLong)
    # Go through all the point pairs and translate them to pixel pairings
    latLonPairs = []
    for point in pixelPairs:
        # Translate the pixel pairs into untranslated points
        ulon = point[0] * gt[1] + gt[0]
        ulat = point[1] * gt[5] + gt[3]
        # Transform the points to the space
        (lon, lat, holder) = ct.TransformPoint(ulon, ulat)
        # Add the point to our return array
        latLonPairs.append([lat, lon])

    return latLonPairs


def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None):
    """The following method translates given latitude/longitude pairs into pixel locations on a given GEOTIF
    This method does not take into account pixel size and assumes a high enough
244
245
246
    image resolution for pixel size to be insignificant

    :param latLonPairs:     The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]]
Daniel Scheffler's avatar
Daniel Scheffler committed
247
248
249
250
251
    :param path_im:         The file location of the input image
    :param geotransform:    GDAL GeoTransform
    :param projection:      GDAL Projection
    :return:            The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]]
    """
252
    assert path_im or geotransform and projection, \
Daniel Scheffler's avatar
Daniel Scheffler committed
253
254
        "GEOP.latLonToPixel: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
    if geotransform and projection:
255
        gt, proj = geotransform, projection
Daniel Scheffler's avatar
Daniel Scheffler committed
256
    else:
257
        ds = gdal.Open(path_im)
Daniel Scheffler's avatar
Daniel Scheffler committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
        gt, proj = ds.GetGeoTransform(), ds.GetProjection()
    # Load the image dataset
    # Create a spatial reference object for the dataset
    srs = osr.SpatialReference()
    srs.ImportFromWkt(proj)
    # Set up the coordinate transformation object
    srsLatLong = srs.CloneGeogCS()
    ct = osr.CoordinateTransformation(srsLatLong, srs)
    # Go through all the point pairs and translate them to latitude/longitude pairings
    pixelPairs = []
    for point in latLonPairs:
        # Change the point locations into the GeoTransform space
        (point[1], point[0], holder) = ct.TransformPoint(point[1], point[0])
        # Translate the x and y coordinates into pixel values
        x = (point[1] - gt[0]) / gt[1]
        y = (point[0] - gt[3]) / gt[5]
        # Add the point to our return array
        pixelPairs.append([int(x), int(y)])
Daniel Scheffler's avatar
Daniel Scheffler committed
276
277
    return pixelPairs

278

279
280
281
282
283
284
285
286
287
288
289
def lonlat_to_pixel(lon, lat, inverse_geo_transform):
    """Translates the given lon, lat to the grid pixel coordinates in data array (zero start)"""
    # transform to utm
    utm_x, utm_y, utm_z = transform_wgs84_to_utm(lon, lat)
    # apply inverse geo tranform
    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
    return pixel_x, abs(pixel_y)  # y pixel is likly a negative value given geo_transform


Daniel Scheffler's avatar
Daniel Scheffler committed
290
291
292
293
294
295
296
297
298
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:
    """
299
300
    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
301
302
303
            for GCP in gcpList]


304
305
306
307
308
309
310
311
312
def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
    """Reprojects any shapely geometry from one projection to another.

    :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>)
    """

    project = partial(pyproj.transform, pyproj.Proj(get_proj4info(proj=prj_src)),
313
                      pyproj.Proj(get_proj4info(proj=prj_tgt)))
314
315
316
317

    return transform(project, shapelyGeometry)  # apply projection


318
def reproject_shapelyPoly(shapelyPoly, tgt_prj):
319
    # type: (Polygon, str) -> Polygon
320
321
322
323
    """Repojects a shapely polygon that has LonLat coordinates into the given target projection.
    :param shapelyPoly:     <shapely.Polygon> the shapely polygon to be reprojected (must have lonlat coordinates)
    :param tgt_prj:         <str> WKT projection string (like GDAL projection)
    """
324
    warnings.warn('reproject_shapelyPoly() is deprecated. Use reproject_shapelyGeometry() instead.', DeprecationWarning)
325
    # TODO nicht sicher, ob wirklich nur lonlat funktioniert
326
327
328

    def get_coordsArr(shpPoly): return np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
    coordsArr = get_coordsArr(shapelyPoly)
329
330
331
    coordsArrTgtPrj = np.zeros_like(coordsArr)

    proj4dict = proj4_to_dict(get_proj4info(proj=tgt_prj))
332
    ReProj = pyproj.Proj(proj=proj4dict['proj'], zone=proj4dict['zone'], ellps=proj4dict['datum'])
333
334
335
336

    x, y = ReProj(coordsArr[:, 0], coordsArr[:, 1], inverse=False)
    coordsArrTgtPrj[:, 0] = x
    coordsArrTgtPrj[:, 1] = y
337
    return Polygon(coordsArrTgtPrj)