Coverage for src/rok4/utils.py: 97%

1"""Provide functions to manipulate OGR / OSR entities 

2""" 

3 

4# -- IMPORTS -- 

5 

6# standard library 

7import os 

8import re 

9from typing import Tuple 

10 

11# 3rd party 

12from osgeo import gdal, ogr, osr 

13 

14# package 

15from rok4.enums import ColorFormat 

16 

17# -- GLOBALS -- 

18ogr.UseExceptions() 

19osr.UseExceptions() 

20gdal.UseExceptions() 

21 

22__SR_BOOK = {} 

23 

24 

25def srs_to_spatialreference(srs: str) -> "osgeo.osr.SpatialReference": 

26 """Convert coordinates system as string to OSR spatial reference 

27 

28 Using a cache, to instanciate a Spatial Reference from a string only once. 

29 

30 Args: 

31 srs (str): coordinates system PROJ4 compliant, with authority and code, like EPSG:3857 or IGNF:LAMB93 

32 

33 Raises: 

34 RuntimeError: Provided SRS is invalid for OSR 

35 

36 Returns: 

37 osgeo.osr.SpatialReference: Corresponding OSR spatial reference 

38 """ 

39 

40 global __SR_BOOK 

41 

42 if srs.upper() not in __SR_BOOK: 

43 authority, code = srs.split(":", 1) 

44 

45 sr = osr.SpatialReference() 

46 if authority.upper() == "EPSG": 

47 sr.ImportFromEPSG(int(code)) 

48 else: 

49 sr.ImportFromProj4(f"+init={srs.upper()} +wktext") 

50 

51 __SR_BOOK[srs.upper()] = sr 

52 

53 return __SR_BOOK[srs.upper()] 

54 

55 

56def bbox_to_geometry( 

57 bbox: Tuple[float, float, float, float], densification: int = 0 

58) -> "osgeo.ogr.Geometry": 

59 """Convert bbox coordinates to OGR geometry 

60 

61 Args: 

62 bbox (Tuple[float, float, float, float]): bounding box (xmin, ymin, xmax, ymax) 

63 densification (int, optional): Number of point to add for each side of bounding box. Defaults to 0. 

64 

65 Raises: 

66 RuntimeError: Provided SRS is invalid for OSR 

67 

68 Returns: 

69 osgeo.ogr.Geometry: Corresponding OGR geometry, with spatial reference if provided 

70 """ 

71 

72 ring = ogr.Geometry(ogr.wkbLinearRing) 

73 

74 if densification > 0: 

75 step_x = (bbox[2] - bbox[0]) / (densification + 1) 

76 step_y = (bbox[3] - bbox[1]) / (densification + 1) 

77 

78 for i in range(densification + 1): 

79 ring.AddPoint(bbox[0] + step_x * i, bbox[1]) 

80 for i in range(densification + 1): 

81 ring.AddPoint(bbox[2], bbox[1] + step_y * i) 

82 for i in range(densification + 1): 

83 ring.AddPoint(bbox[2] - step_x * i, bbox[3]) 

84 for i in range(densification + 1): 

85 ring.AddPoint(bbox[0], bbox[3] - step_y * i) 

86 ring.AddPoint(bbox[0], bbox[1]) 

87 

88 else: 

89 ring.AddPoint(bbox[0], bbox[1]) 

90 ring.AddPoint(bbox[2], bbox[1]) 

91 ring.AddPoint(bbox[2], bbox[3]) 

92 ring.AddPoint(bbox[0], bbox[3]) 

93 ring.AddPoint(bbox[0], bbox[1]) 

94 

95 geom = ogr.Geometry(ogr.wkbPolygon) 

96 geom.AddGeometry(ring) 

97 geom.SetCoordinateDimension(2) 

98 

99 return geom 

100 

101 

102def reproject_bbox( 

103 bbox: Tuple[float, float, float, float], srs_src: str, srs_dst: str, densification: int = 5 

104) -> Tuple[float, float, float, float]: 

105 """Return bounding box in other coordinates system 

106 

107 Points are added to be sure output bounding box contains input bounding box 

108 

109 Args: 

110 bbox (Tuple[float, float, float, float]): bounding box (xmin, ymin, xmax, ymax) with source coordinates system 

111 srs_src (str): source coordinates system 

112 srs_dst (str): destination coordinates system 

113 densification (int, optional): Number of point to add for each side of bounding box. Defaults to 5. 

114 

115 Returns: 

116 Tuple[float, float, float, float]: bounding box (xmin, ymin, xmax, ymax) with destination coordinates system 

117 """ 

118 

119 sr_src = srs_to_spatialreference(srs_src) 

120 sr_src_inv = sr_src.EPSGTreatsAsLatLong() or sr_src.EPSGTreatsAsNorthingEasting() 

121 

122 sr_dst = srs_to_spatialreference(srs_dst) 

123 sr_dst_inv = sr_dst.EPSGTreatsAsLatLong() or sr_dst.EPSGTreatsAsNorthingEasting() 

124 

125 if sr_src.IsSame(sr_dst) and sr_src_inv == sr_dst_inv: 

126 # Les système sont vraiment les même, avec le même ordre des axes 

127 return bbox 

128 elif sr_src.IsSame(sr_dst) and sr_src_inv != sr_dst_inv: 

129 # Les système sont les même pour OSR, mais l'ordre des axes est différent 

130 return (bbox[1], bbox[0], bbox[3], bbox[2]) 

131 

132 # Systèmes différents 

133 

134 bbox_src = bbox_to_geometry(bbox, densification) 

135 bbox_src.AssignSpatialReference(sr_src) 

136 

137 bbox_dst = bbox_src.Clone() 

138 os.environ["OGR_ENABLE_PARTIAL_REPROJECTION"] = "YES" 

139 bbox_dst.TransformTo(sr_dst) 

140 

141 env = bbox_dst.GetEnvelope() 

142 return (env[0], env[2], env[1], env[3]) 

143 

144 

145def reproject_point( 

146 point: Tuple[float, float], 

147 sr_src: "osgeo.osr.SpatialReference", 

148 sr_dst: "osgeo.osr.SpatialReference", 

149) -> Tuple[float, float]: 

150 """Reproject a point 

151 

152 Args: 

153 point (Tuple[float, float]): source spatial reference point 

154 sr_src (osgeo.osr.SpatialReference): source spatial reference 

155 sr_dst (osgeo.osr.SpatialReference): destination spatial reference 

156 

157 Returns: 

158 Tuple[float, float]: X/Y in destination spatial reference 

159 """ 

160 

161 sr_src_inv = sr_src.EPSGTreatsAsLatLong() or sr_src.EPSGTreatsAsNorthingEasting() 

162 sr_dst_inv = sr_dst.EPSGTreatsAsLatLong() or sr_dst.EPSGTreatsAsNorthingEasting() 

163 

164 if sr_src.IsSame(sr_dst) and sr_src_inv == sr_dst_inv: 

165 # Les système sont vraiment les mêmes, avec le même ordre des axes 

166 return (point[0], point[1]) 

167 elif sr_src.IsSame(sr_dst) and sr_src_inv != sr_dst_inv: 

168 # Les système sont les même pour OSR, mais l'ordre des axes est différent 

169 return (point[1], point[0]) 

170 

171 # Systèmes différents 

172 ct = osr.CreateCoordinateTransformation(sr_src, sr_dst) 

173 x_dst, y_dst, z_dst = ct.TransformPoint(point[0], point[1]) 

174 

175 return (x_dst, y_dst) 

176 

177 

178def compute_bbox(source_dataset: gdal.Dataset) -> Tuple: 

179 """Image boundingbox computing method 

180 

181 Args: 

182 source_dataset (gdal.Dataset): Dataset instanciated 

183 from the raster image 

184 

185 Limitations: 

186 Image's axis must be parallel to SRS' axis 

187 

188 Raises: 

189 AttributeError: source_dataset is not a gdal.Dataset instance. 

190 Exception: The dataset does not contain transform data. 

191 """ 

192 bbox = None 

193 transform_vector = source_dataset.GetGeoTransform() 

194 if transform_vector is None: 

195 raise Exception( 

196 "No transform vector found in the dataset created from " 

197 + f"the following file : {source_dataset.GetFileList()[0]}" 

198 ) 

199 width = source_dataset.RasterXSize 

200 height = source_dataset.RasterYSize 

201 x_range = ( 

202 transform_vector[0], 

203 transform_vector[0] + width * transform_vector[1] + height * transform_vector[2], 

204 ) 

205 y_range = ( 

206 transform_vector[3], 

207 transform_vector[3] + width * transform_vector[4] + height * transform_vector[5], 

208 ) 

209 spatial_ref = source_dataset.GetSpatialRef() 

210 if spatial_ref is not None and spatial_ref.GetDataAxisToSRSAxisMapping() == [2, 1]: 

211 # Coordonnées terrain de type (latitude, longitude) 

212 # => on permute les coordonnées terrain par rapport à l'image 

213 bbox = (min(y_range), min(x_range), max(y_range), max(x_range)) 

214 else: 

215 # Coordonnées terrain de type (longitude, latitude) ou pas de SRS 

216 # => les coordonnées terrain sont dans le même ordre que celle de l'image 

217 bbox = (min(x_range), min(y_range), max(x_range), max(y_range)) 

218 return bbox 

219 

220 

221def compute_format(dataset: gdal.Dataset, path: str = None) -> ColorFormat: 

222 """Image color format computing method 

223 

224 Args: 

225 dataset (gdal.Dataset): Dataset instanciated from the image 

226 path (str, optionnal): path to the original file/object 

227 

228 Raises: 

229 AttributeError: source_dataset is not a gdal.Dataset instance. 

230 Exception: No color band found or unsupported color format. 

231 """ 

232 color_format = None 

233 if path is None: 

234 path = dataset.GetFileList()[0] 

235 if dataset.RasterCount < 1: 

236 raise Exception(f"Image {path} contains no color band.") 

237 

238 band_1_datatype = dataset.GetRasterBand(1).DataType 

239 data_type_name = gdal.GetDataTypeName(band_1_datatype) 

240 data_type_size = gdal.GetDataTypeSize(band_1_datatype) 

241 color_interpretation = dataset.GetRasterBand(1).GetRasterColorInterpretation() 

242 color_name = None 

243 if color_interpretation is not None: 

244 color_name = gdal.GetColorInterpretationName(color_interpretation) 

245 compression_regex_match = re.search(r"COMPRESSION\s*=\s*PACKBITS", gdal.Info(dataset)) 

246 

247 if ( 

248 data_type_name == "Byte" 

249 and data_type_size == 8 

250 and color_name == "Palette" 

251 and compression_regex_match 

252 ): 

253 # Compris par libTIFF comme du noir et blanc sur 1 bit 

254 color_format = ColorFormat.BIT 

255 elif data_type_name == "Byte" and data_type_size == 8: 

256 color_format = ColorFormat.UINT8 

257 elif data_type_name == "Float32" and data_type_size == 32: 

258 color_format = ColorFormat.FLOAT32 

259 else: 

260 raise Exception( 

261 f"Unsupported color format for image {path} : " 

262 + f"'{data_type_name}' ({data_type_size} bits)" 

263 ) 

264 return color_format