Coverage for py_tools_ds/geo/coord_trafo.py: 78%

1# -*- coding: utf-8 -*- 

2 

3# py_tools_ds - A collection of geospatial data analysis tools that simplify standard 

4# operations when handling geospatial raster and vector data as well as projections. 

5# 

6# Copyright (C) 2016-2021 

7# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

8# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam, 

9# Germany (https://www.gfz-potsdam.de/) 

10# 

11# This software was developed within the context of the GeoMultiSens project funded 

12# by the German Federal Ministry of Education and Research 

13# (project grant code: 01 IS 14 010 A-C). 

14# 

15# Licensed under the Apache License, Version 2.0 (the "License"); 

16# you may not use this file except in compliance with the License. 

17# You may obtain a copy of the License at 

18# 

19# http://www.apache.org/licenses/LICENSE-2.0 

20# 

21# Unless required by applicable law or agreed to in writing, software 

22# distributed under the License is distributed on an "AS IS" BASIS, 

23# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

24# See the License for the specific language governing permissions and 

25# limitations under the License. 

26 

27import pyproj 

28import numpy as np 

29from shapely.ops import transform 

30from typing import Union, Sequence, List # noqa F401 # flake8 issue 

31from osgeo import gdal, osr 

32 

33__author__ = "Daniel Scheffler" 

34 

35 

36def transform_utm_to_wgs84(easting, northing, zone, south=False): 

37 # type: (float, float, int, bool) -> (float, float) 

38 """Transform an UTM coordinate to lon, lat, altitude.""" 

39 UTM = pyproj.Proj(proj='utm', 

40 zone=abs(zone), 

41 ellps='WGS84', 

42 south=(zone < 0 or south)) 

43 return UTM(easting, northing, inverse=True) 

44 

45 

46def get_utm_zone(longitude): 

47 # type: (float) -> int 

48 """Get the UTM zone corresponding to the given longitude.""" 

49 return int(1 + (longitude + 180.0) / 6.0) 

50 

51 

52def transform_wgs84_to_utm(lon, lat, zone=None): 

53 # type: (float, float, int) -> (float, float) 

54 """Return easting, northing, altitude for the given Lon/Lat coordinate. 

55 

56 :param lon: longitude 

57 :param lat: latitude 

58 :param zone: UTM zone (if not set it is derived automatically) 

59 """ 

60 if not (-180. <= lon <= 180. and -90. <= lat <= 90.): 

61 raise ValueError((lon, lat), 'Input coordinates for transform_wgs84_to_utm() out of range.') 

62 

63 UTM = pyproj.Proj(proj='utm', 

64 zone=get_utm_zone(lon) if zone is None else zone, 

65 ellps='WGS84') 

66 return UTM(lon, lat) 

67 

68 

69def transform_any_prj(prj_src, prj_tgt, x, y): 

70 # type: (Union[str, int], Union[str, int], Union[float, np.ndarray], Union[float, np.ndarray]) -> (Union[float, np.ndarray], Union[float, np.ndarray]) # noqa 

71 """Transform X/Y data from any source projection to any target projection. 

72 

73 :param prj_src: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>) 

74 :param prj_tgt: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>) 

75 :param x: X-coordinate(s) to be transformed (either <float> or 1D-<np.ndarray>) 

76 :param y: X-coordinate(s) to be transformed (either <float> or 1D-<np.ndarray>) 

77 :return: 

78 """ 

79 transformer = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(prj_src), 

80 pyproj.CRS.from_user_input(prj_tgt), 

81 always_xy=True) 

82 return transformer.transform(x, y) 

83 

84 

85def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None): 

86 # type: (str, str, np.ndarray, np.ndarray, np.ndarray) -> Sequence[np.ndarray] 

87 """Transform a geolocation array from one projection into another. 

88 

89 HINT: This function is faster than transform_any_prj but works only for geolocation arrays. 

90 

91 :param prj_src: WKT string 

92 :param prj_tgt: WKT string 

93 :param Xarr: np.ndarray of shape (rows,cols) 

94 :param Yarr: np.ndarray of shape (rows,cols) 

95 :param Zarr: np.ndarray of shape (rows,cols) 

96 :return: 

97 """ 

98 drv = gdal.GetDriverByName('MEM') 

99 geoloc_ds = drv.Create('geoloc', Xarr.shape[1], Xarr.shape[0], 3, gdal.GDT_Float64) 

100 geoloc_ds.GetRasterBand(1).WriteArray(Xarr) 

101 geoloc_ds.GetRasterBand(2).WriteArray(Yarr) 

102 if Zarr is not None: 

103 geoloc_ds.GetRasterBand(3).WriteArray(Zarr) 

104 

105 def strip_wkt(wkt): 

106 return wkt\ 

107 .replace('\n', '')\ 

108 .replace('\r', '')\ 

109 .replace(' ', '') 

110 

111 transformer = gdal.Transformer(None, None, ['SRC_SRS=' + strip_wkt(prj_src), 

112 'DST_SRS=' + strip_wkt(prj_tgt)]) 

113 status = transformer.TransformGeolocations(geoloc_ds.GetRasterBand(1), 

114 geoloc_ds.GetRasterBand(2), 

115 geoloc_ds.GetRasterBand(3)) 

116 

117 if status: 

118 raise Exception('Error transforming coordinate array: ' + gdal.GetLastErrorMsg()) 

119 

120 Xarr = geoloc_ds.GetRasterBand(1).ReadAsArray() 

121 Yarr = geoloc_ds.GetRasterBand(2).ReadAsArray() 

122 

123 if Zarr is not None: 

124 Zarr = geoloc_ds.GetRasterBand(3).ReadAsArray() 

125 return Xarr, Yarr, Zarr 

126 else: 

127 return Xarr, Yarr 

128 

129 

130def mapXY2imXY(mapXY, gt): 

131 # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray] 

132 """Translate the given geo coordinates into pixel locations according to the given image geotransform. 

133 

134 :param mapXY: <tuple, np.ndarray> The geo coordinates to be translated in the form (x,y) or as np.ndarray [Nx1]. 

135 :param gt: <list> GDAL geotransform 

136 :returns: <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2] 

137 """ 

138 if isinstance(mapXY, np.ndarray): 

139 ndim = mapXY.ndim 

140 if ndim == 1: 

141 mapXY = mapXY.reshape(1, 2) 

142 assert mapXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % mapXY.shape 

143 imXY = np.empty_like(mapXY, dtype=float) 

144 imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1] 

145 imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5] 

146 return imXY if ndim > 1 else imXY.flatten() 

147 else: 

148 return (mapXY[0] - gt[0]) / gt[1],\ 

149 (mapXY[1] - gt[3]) / gt[5] 

150 

151 

152def imXY2mapXY(imXY, gt): 

153 # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray] 

154 """Translate the given pixel X/Y locations into geo coordinates according to the given image geotransform. 

155 

156 :param imXY: <tuple, np.ndarray> The image coordinates to be translated in the form (x,y) or as np.ndarray [Nx1]. 

157 :param gt: <list> GDAL geotransform 

158 :returns: <tuple, np.ndarray> geo coordinate tuple X/Y (mapX, mapY) or np.ndarray [Nx2] 

159 """ 

160 if isinstance(imXY, np.ndarray): 

161 ndim = imXY.ndim 

162 if imXY.ndim == 1: 

163 imXY = imXY.reshape(1, 2) 

164 assert imXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % imXY.shape 

165 mapXY = np.empty_like(imXY, dtype=float) 

166 mapXY[:, 0] = gt[0] + imXY[:, 0] * abs(gt[1]) 

167 mapXY[:, 1] = gt[3] - imXY[:, 1] * abs(gt[5]) 

168 return mapXY if ndim > 1 else mapXY.flatten() 

169 else: 

170 return (gt[0] + imXY[0] * abs(gt[1])),\ 

171 (gt[3] - imXY[1] * abs(gt[5])) 

172 

173 

174def mapYX2imYX(mapYX, gt): 

175 return (mapYX[0] - gt[3]) / gt[5],\ 

176 (mapYX[1] - gt[0]) / gt[1] 

177 

178 

179def imYX2mapYX(imYX, gt): 

180 return gt[3] - (imYX[0] * abs(gt[5])),\ 

181 gt[0] + (imYX[1] * gt[1]) 

182 

183 

184def pixelToMapYX(pixelCoords, geotransform, projection): 

185 # MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform) 

186 

187 # Create a spatial reference object 

188 srs = osr.SpatialReference() 

189 srs.ImportFromWkt(projection) 

190 

191 # Set up the coordinate transformation object 

192 ct = osr.CoordinateTransformation(srs, srs) 

193 

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

195 

196 mapXYPairs = [ct.TransformPoint(pixX * geotransform[1] + geotransform[0], 

197 pixY * geotransform[5] + geotransform[3])[:2] 

198 for pixX, pixY in pixelCoords] 

199 

200 return [[mapY, mapX] for mapX, mapY in mapXYPairs] 

201 

202 

203def pixelToLatLon(pixelPairs, geotransform, projection): 

204 # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]] 

205 """Translate the given pixel X/Y locations into latitude/longitude locations. 

206 

207 :param pixelPairs: Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]] 

208 :param geotransform: GDAL GeoTransform 

209 :param projection: WKT string 

210 :returns: The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]] 

211 """ 

212 MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform) 

213 lons, lats = transform_any_prj(projection, 4326, MapXYPairs[:, 0], MapXYPairs[:, 1]) 

214 

215 # validate output lat/lon values 

216 if not (-180 <= np.min(lons) <= 180) or not (-180 <= np.max(lons) <= 180): 

217 raise RuntimeError('Output longitude values out of bounds.') 

218 if not (-90 <= np.min(lats) <= 90) or not (-90 <= np.max(lats) <= 90): 

219 raise RuntimeError('Output latitudes values out of bounds.') 

220 

221 LatLonPairs = np.vstack([lats, lons]).T.tolist() 

222 

223 return LatLonPairs 

224 

225 

226def latLonToPixel(latLonPairs, geotransform, projection): 

227 # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]] 

228 """Translate the given latitude/longitude pairs into pixel X/Y locations. 

229 

230 :param latLonPairs: The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]] 

231 :param geotransform: GDAL GeoTransform 

232 :param projection: WKT string 

233 :return: The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]] 

234 """ 

235 latLonArr = np.array(latLonPairs) 

236 

237 # validate input lat/lon values 

238 if not (-180 <= np.min(latLonArr[:, 1]) <= 180) or not (-180 <= np.max(latLonArr[:, 1]) <= 180): 

239 raise ValueError('Longitude value out of bounds.') 

240 if not (-90 <= np.min(latLonArr[:, 0]) <= 90) or not (-90 <= np.max(latLonArr[:, 0]) <= 90): 

241 raise ValueError('Latitude value out of bounds.') 

242 

243 mapXs, mapYs = transform_any_prj(4326, projection, latLonArr[:, 1], latLonArr[:, 0]) 

244 imXYarr = mapXY2imXY(np.vstack([mapXs, mapYs]).T, geotransform) 

245 

246 pixelPairs = imXYarr.tolist() 

247 

248 return pixelPairs 

249 

250 

251def lonlat_to_pixel(lon, lat, inverse_geo_transform): 

252 """Translate the given lon, lat to the grid pixel coordinates in data array (zero start).""" 

253 # transform to utm 

254 utm_x, utm_y = transform_wgs84_to_utm(lon, lat) 

255 

256 # apply inverse geo transform 

257 pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y) 

258 pixel_x = int(pixel_x) - 1 # adjust to 0 start for array 

259 pixel_y = int(pixel_y) - 1 # adjust to 0 start for array 

260 

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

262 

263 

264def transform_GCPlist(gcpList, prj_src, prj_tgt): 

265 """ 

266 

267 :param gcpList: <list> list of ground control points in the output coordinate system 

268 to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...] 

269 :param prj_src: WKT string 

270 :param prj_tgt: WKT string 

271 :return: 

272 """ 

273 return [gdal.GCP(*(list(transform_any_prj(prj_src, prj_tgt, GCP.GCPX, GCP.GCPY)) + 

274 [GCP.GCPZ, GCP.GCPPixel, GCP.GCPLine])) # Python 2.7 compatible syntax 

275 for GCP in gcpList] 

276 

277 

278def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt): 

279 """Reproject any shapely geometry from one projection to another. 

280 

281 :param shapelyGeometry: any shapely geometry instance 

282 :param prj_src: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>) 

283 :param prj_tgt: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>) 

284 """ 

285 project = pyproj.Transformer.from_proj( 

286 pyproj.CRS.from_user_input(prj_src), 

287 pyproj.CRS.from_user_input(prj_tgt), 

288 always_xy=True) 

289 

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