Coverage for py_tools_ds/geo/vector/conversion.py: 47%

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 numpy as np 

28 

29# custom 

30from shapely.geometry import shape, mapping, box 

31from shapely.geometry import Polygon # noqa F401 # flake8 issue 

32from shapely import wkt 

33from pyproj import CRS 

34from osgeo import gdal, ogr, osr 

35 

36from ..coord_trafo import imYX2mapYX, mapYX2imYX, pixelToMapYX 

37from ...dtypes.conversion import get_dtypeStr, dTypeDic_NumPy2GDAL 

38 

39__author__ = "Daniel Scheffler" 

40 

41 

42def shapelyImPoly_to_shapelyMapPoly_withPRJ(shapelyImPoly, gt, prj): 

43 # ACTUALLY PRJ IS NOT NEEDED BUT THIS FUNCTION RETURNS OTHER VALUES THAN shapelyImPoly_to_shapelyMapPoly 

44 geojson = mapping(shapelyImPoly) 

45 coords = list(geojson['coordinates'][0]) 

46 coordsYX = pixelToMapYX(coords, geotransform=gt, projection=prj) 

47 coordsXY = tuple([(i[1], i[0]) for i in coordsYX]) 

48 geojson['coordinates'] = (coordsXY,) 

49 return shape(geojson) 

50 

51 

52def shapelyImPoly_to_shapelyMapPoly(shapelyBox, gt): 

53 xmin, ymin, xmax, ymax = shapelyBox.bounds 

54 ymax, xmin = imYX2mapYX((ymax, xmin), gt) 

55 ymin, xmax = imYX2mapYX((ymin, xmax), gt) 

56 return box(xmin, ymin, xmax, ymax) 

57 

58 

59def shapelyBox2BoxYX(shapelyBox, coord_type='image'): 

60 xmin, ymin, xmax, ymax = shapelyBox.bounds 

61 assert coord_type in ['image', 'map'] 

62 

63 if coord_type == 'image': 

64 UL_YX = ymin, xmin 

65 UR_YX = ymin, xmax 

66 LR_YX = ymax, xmax 

67 LL_YX = ymax, xmin 

68 else: 

69 UL_YX = ymax, xmin 

70 UR_YX = ymax, xmax 

71 LR_YX = ymin, xmax 

72 LL_YX = ymin, xmin 

73 

74 return UL_YX, UR_YX, LR_YX, LL_YX 

75 

76 

77def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt): 

78 # type: (Polygon, tuple) -> list 

79 """Convert each vertex coordinate of a shapely polygon into image coordinates corresponding to the given 

80 geotransform without respect to invalid image coordinates. Those must be filtered later. 

81 

82 :param shapelyPoly: <shapely.Polygon> 

83 :param im_gt: <list> the GDAL geotransform of the target image 

84 """ 

85 def get_coordsArr(shpPoly): return np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1) 

86 coordsArr = get_coordsArr(shapelyPoly) 

87 boxImXY = [mapYX2imYX((Y, X), im_gt) for X, Y in coordsArr.tolist()] # FIXME incompatible to GMS version 

88 boxImXY = [(i[1], i[0]) for i in boxImXY] 

89 return boxImXY 

90 

91 

92def round_shapelyPoly_coords(shapelyPoly, precision=10): 

93 # type: (Polygon, int) -> Polygon 

94 """Round the coordinates of the given shapely polygon. 

95 

96 :param shapelyPoly: the shapely polygon 

97 :param precision: number of decimals 

98 :return: 

99 """ 

100 return wkt.loads(wkt.dumps(shapelyPoly, rounding_precision=precision)) 

101 

102 

103def points_to_raster(points, values, tgt_res, prj=None, fillVal=None): 

104 # type: (np.ndarray, np.ndarray, float, str, float) -> (np.ndarray, list, str) 

105 """ 

106 Convert a set of point geometries with associated values into a raster array. 

107 

108 :param points: list or 1D numpy.ndarray containings shapely.geometry point geometries 

109 :param values: list or 1D numpy.ndarray containing int or float values 

110 :param tgt_res: target resolution in projection units 

111 :param prj: WKT projection string 

112 :param fillVal: fill value used to fill in where no point geometry is available 

113 """ 

114 ds = ogr.GetDriverByName("Memory").CreateDataSource('wrk') 

115 

116 if prj is not None: 

117 if int(gdal.__version__[0]) < 3: 

118 # noinspection PyTypeChecker 

119 prj = CRS(prj).to_wkt(version="WKT1_GDAL") 

120 

121 srs = osr.SpatialReference() 

122 srs.ImportFromWkt(prj) 

123 else: 

124 srs = None 

125 

126 layer = ds.CreateLayer('', srs, ogr.wkbPoint) 

127 

128 # create field 

129 DTypeStr = get_dtypeStr(values if isinstance(values, np.ndarray) else values[0]) 

130 FieldType = ogr.OFTInteger if DTypeStr.startswith('int') else ogr.OFTReal 

131 FieldDefn = ogr.FieldDefn('VAL', FieldType) 

132 if DTypeStr.startswith('float'): 

133 FieldDefn.SetPrecision(6) 

134 layer.CreateField(FieldDefn) # Add one attribute 

135 

136 for i in range(len(points)): 

137 # Create a new feature (attribute and geometry) 

138 feat = ogr.Feature(layer.GetLayerDefn()) 

139 feat.SetGeometry(ogr.CreateGeometryFromWkb(points[i].wkb)) # Make a geometry, from Shapely object 

140 feat.SetField('VAL', values[i]) 

141 

142 layer.CreateFeature(feat) 

143 feat.Destroy() 

144 

145 x_min, x_max, y_min, y_max = layer.GetExtent() 

146 

147 # Create the destination data source 

148 cols = int((x_max - x_min) / tgt_res) 

149 rows = int((y_max - y_min) / tgt_res) 

150 target_ds = gdal.GetDriverByName('MEM').Create('raster', cols, rows, 1, dTypeDic_NumPy2GDAL[DTypeStr]) 

151 target_ds.SetGeoTransform((x_min, tgt_res, 0, y_max, 0, -tgt_res)) 

152 target_ds.SetProjection(prj if prj else '') 

153 band = target_ds.GetRasterBand(1) 

154 if fillVal is not None: 

155 band.Fill(fillVal) 

156 band.FlushCache() 

157 

158 # Rasterize 

159 gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=VAL"]) 

160 

161 out_arr = target_ds.GetRasterBand(1).ReadAsArray() 

162 out_gt = target_ds.GetGeoTransform() 

163 out_prj = target_ds.GetProjection() 

164 

165 del target_ds, ds, layer, band 

166 

167 return out_arr, out_gt, out_prj