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
from shapely.ops import transform
Daniel Scheffler's avatar
Daniel Scheffler committed
10
from typing import Union  # noqa F401  # flake8 issue
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
__author__ = "Daniel Scheffler"

Daniel Scheffler's avatar
Daniel Scheffler committed
23 24 25 26 27 28 29 30 31 32 33

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
34
    #    utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system)
Daniel Scheffler's avatar
Daniel Scheffler committed
35 36
    #    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]
37
    UTM = pyproj.Proj(proj='utm', zone=abs(zone), ellps='WGS84', south=(zone < 0 or south))
Daniel Scheffler's avatar
Daniel Scheffler committed
38 39 40 41 42 43 44 45 46
    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:
    """
47 48 49 50
    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
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
    #    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
67
    #    wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system)
Daniel Scheffler's avatar
Daniel Scheffler committed
68 69 70 71 72
    #    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)


73
def transform_any_prj(prj_src, prj_tgt, x, y):
74 75 76 77 78 79 80 81
    """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
82 83
    prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
    prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
84 85
    x, y = pyproj.transform(prj_src, prj_tgt, x, y)
    return x, y
Daniel Scheffler's avatar
Daniel Scheffler committed
86 87


88 89 90 91 92 93 94 95 96 97 98
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:
    """
99 100
    drv = gdal.GetDriverByName('MEM')
    geoloc_ds = drv.Create('geoloc', Xarr.shape[1], Xarr.shape[0], 3, gdal.GDT_Float64)
101 102 103 104 105 106
    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])
107
    status = transformer.TransformGeolocations(geoloc_ds.GetRasterBand(1), geoloc_ds.GetRasterBand(2),
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
                                               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
123
def mapXY2imXY(mapXY, gt):
Daniel Scheffler's avatar
Daniel Scheffler committed
124
    # type: (tuple, Union[list, tuple]) -> Union[tuple, np.ndarray]
125
    """Translates given geo coordinates into pixel locations according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
126

127
    :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
128
    :param gt:      <list> GDAL geotransform
129
    :returns:       <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2]
Daniel Scheffler's avatar
Daniel Scheffler committed
130
    """
131
    if isinstance(mapXY, np.ndarray):
132 133 134
        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
135
        imXY = np.empty_like(mapXY)
136 137
        imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1]
        imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5]
138 139
        return imXY
    else:
140
        return (mapXY[0] - gt[0]) / gt[1], (mapXY[1] - gt[3]) / gt[5]
Daniel Scheffler's avatar
Daniel Scheffler committed
141 142 143


def imXY2mapXY(imXY, gt):
Daniel Scheffler's avatar
Daniel Scheffler committed
144
    # type: (tuple, Union[list, tuple]) -> Union[tuple, np.ndarray]
145
    """Translates given pixel locations into geo coordinates according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
146

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


163 164 165 166 167 168
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])
169 170 171 172 173 174


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:
175
        gt, proj = geotransform, projection
176
    else:
177
        ds = gdal.Open(path_im)
178 179 180 181
        gt, proj = ds.GetGeoTransform(), ds.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(proj)
    # Set up the coordinate transformation object
182
    ct = osr.CoordinateTransformation(srs, srs)
183 184

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

    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
199 200 201 202 203 204 205 206 207 208 209 210 211 212
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:
213
        gt, proj = geotransform, projection
Daniel Scheffler's avatar
Daniel Scheffler committed
214
    else:
215
        ds = gdal.Open(path_im)
Daniel Scheffler's avatar
Daniel Scheffler committed
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
        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
241 242 243
    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
244 245 246 247 248
    :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]]
    """
249
    assert path_im or geotransform and projection, \
Daniel Scheffler's avatar
Daniel Scheffler committed
250 251
        "GEOP.latLonToPixel: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
    if geotransform and projection:
252
        gt, proj = geotransform, projection
Daniel Scheffler's avatar
Daniel Scheffler committed
253
    else:
254
        ds = gdal.Open(path_im)
Daniel Scheffler's avatar
Daniel Scheffler committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
        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
273 274
    return pixelPairs

275

276 277 278 279 280 281 282 283 284 285 286
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
287 288 289 290 291 292 293 294 295
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:
    """
296 297
    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
298 299 300
            for GCP in gcpList]


301 302 303 304 305 306 307 308 309
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)),
310
                      pyproj.Proj(get_proj4info(proj=prj_tgt)))
311 312 313 314

    return transform(project, shapelyGeometry)  # apply projection


315
def reproject_shapelyPoly(shapelyPoly, tgt_prj):
316
    # type: (Polygon, str) -> Polygon
317 318 319 320
    """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)
    """
321
    warnings.warn('reproject_shapelyPoly() is deprecated. Use reproject_shapelyGeometry() instead.', DeprecationWarning)
322
    # TODO nicht sicher, ob wirklich nur lonlat funktioniert
323 324 325

    def get_coordsArr(shpPoly): return np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
    coordsArr = get_coordsArr(shapelyPoly)
326 327 328
    coordsArrTgtPrj = np.zeros_like(coordsArr)

    proj4dict = proj4_to_dict(get_proj4info(proj=tgt_prj))
329
    ReProj = pyproj.Proj(proj=proj4dict['proj'], zone=proj4dict['zone'], ellps=proj4dict['datum'])
330 331 332 333

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