coord_trafo.py 10.8 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"

import pyproj
5
6
import numpy as np
from shapely.geometry import Polygon
Daniel Scheffler's avatar
Daniel Scheffler committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

try:
    from osgeo import ogr
    from osgeo import osr
    from osgeo import gdal
    from osgeo import gdalnumeric
    from osgeo.gdalconst import *
except ImportError:
    import ogr
    import osr
    import gdal
    import gdalnumeric
    from gdalconst import *


22
from .projection import get_proj4info, proj4_to_dict
Daniel Scheffler's avatar
Daniel Scheffler committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49



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
    #    utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (, )
    #    print('hier2')
    #    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]
    UTM   = pyproj.Proj(proj='utm', zone=abs(zone), ellps='WGS84',south=(zone<0 or south))
    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
    assert -180.<=lon<=180. and -90.<=lat<=90., 'Input coordinates for transform_wgs84_to_utm() out of range. ' \
                                                'Got %s/%s.' %(lon,lat)
Daniel Scheffler's avatar
Daniel Scheffler committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
    get_utm_zone = lambda longitude: int(1 + (longitude + 180.0) / 6.0)
    #    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
    #    wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (, )
    #    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)


def transform_any_prj(prj_src,prj_tgt,x,y):
    prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
    prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
    x,y     = pyproj.transform(prj_src,prj_tgt,x,y)
    return x,y


def mapXY2imXY(mapXY, gt):
    # type: (tuple,list) -> tuple
84
    """Translates given geo coordinates into pixel locations according to the given image geotransform.
Daniel Scheffler's avatar
Daniel Scheffler committed
85

86
    :param mapXY:   <tuple> The geo coordinates to be translated in the form (x,y).
Daniel Scheffler's avatar
Daniel Scheffler committed
87
88
89
90
91
92
93
94
    :param gt:      <list> GDAL geotransform
    :returns:       <tuple> image coordinate tuple X/Y (column, row)
    """
    return (mapXY[0]-gt[0])/gt[1],(mapXY[1]-gt[3])/gt[5]


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

    :param imXY:    <tuple> The image coordinates to be translated in the form (x,y).
    :param gt:      <list> GDAL geotransform
99
    :returns:       <tuple> geo coordinate tuple X/Y (mapX, mapY)
Daniel Scheffler's avatar
Daniel Scheffler committed
100
101
102
103
    """
    return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5]))


104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
mapYX2imYX = lambda mapYX,gt: ((mapYX[0]-gt[3])/gt[5],     (mapYX[1]-gt[0])/gt[1])
imYX2mapYX = lambda imYX ,gt: (gt[3]-(imYX[0]*abs(gt[5])), gt[0]+(imYX[1]*gt[1]))


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:
        gt,proj = geotransform, projection
    else:
        ds       = gdal.Open(path_im)
        gt, proj = ds.GetGeoTransform(), ds.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(proj)
    # Set up the coordinate transformation object
    ct  = osr.CoordinateTransformation(srs, srs)

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

    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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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:
        gt,proj = geotransform, projection
    else:
        ds       = gdal.Open(path_im)
        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
	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]]
    :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]]
    """
    assert path_im or geotransform  and projection, \
        "GEOP.latLonToPixel: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
    if geotransform and projection:
        gt,proj = geotransform, projection
    else:
        ds       = gdal.Open(path_im)
        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
209
210
211
    return pixelPairs


212
213
214
215
216
217
218
219
220
221
222
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
223
224
225
226
227
228
229
230
231
232
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:
    """
    return [gdal.GCP(*transform_any_prj(prj_src, prj_tgt, GCP.GCPX, GCP.GCPY), GCP.GCPZ, GCP.GCPPixel, GCP.GCPLine)
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
            for GCP in gcpList]


def reproject_shapelyPoly(shapelyPoly, tgt_prj):
    # type: (shapely.Polygon, str) -> shapely.Polygon
    """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)
    """
    # TODO nicht sicher, ob wirklich nur lonlat funktioniert
    get_coordsArr   = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
    coordsArr       = get_coordsArr(shapelyPoly)
    coordsArrTgtPrj = np.zeros_like(coordsArr)

    proj4dict = proj4_to_dict(get_proj4info(proj=tgt_prj))
    ReProj    = pyproj.Proj(proj=proj4dict['proj'], zone=proj4dict['zone'], ellps=proj4dict['datum'])

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