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):
Daniel Scheffler's avatar
Daniel Scheffler committed
127
    # type: (tuple, Union[list, tuple]) -> _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):
Daniel Scheffler's avatar
Daniel Scheffler committed
147
    # type: (tuple, Union[list, tuple]) -> _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)