Coverage for gws-app/gws/lib/gdalx/__init__.py: 89%

444 statements  

« prev     ^ index     » next       coverage.py v7.15.4, created at 2026-08-24 12:46 +0200

1"""GDAL/OGR wrapper.""" 

2 

3from typing import Any, Optional, Iterable, cast 

4 

5import datetime 

6import decimal 

7import contextlib 

8import numpy as np 

9 

10from osgeo import gdal 

11from osgeo import ogr 

12from osgeo import osr 

13 

14import gws 

15import gws.base.shape 

16import gws.lib.crs 

17import gws.lib.bounds 

18import gws.lib.image 

19import gws.lib.datetimex as datetimex 

20 

21 

22class Error(gws.Error): 

23 pass 

24 

25 

26class DriverInfo(gws.Data): 

27 index: int 

28 name: str 

29 longName: str 

30 extensions: list[str] 

31 metaData: dict 

32 

33 

34def get_drivers() -> list[DriverInfo]: 

35 """Enumerate GDAL drivers.""" 

36 

37 di = gws.u.get_app_global('gdal_drivers', _fetch_driver_infos) 

38 return di.infos 

39 

40 

41def get_driver(name: str) -> Optional[DriverInfo]: 

42 """Get driver info by name.""" 

43 

44 for di in get_drivers(): 

45 if di.name == name: 

46 return di 

47 

48 

49def supported_attribute_types(): 

50 return list(_ATTR_TO_OGR.keys()) 

51 

52 

53@contextlib.contextmanager 

54def gdal_config(options: dict): 

55 """Temporarily set GDAL config options.""" 

56 

57 prev = {} 

58 for key, value in options.items(): 

59 prev[key] = gdal.GetConfigOption(key) 

60 gdal.SetConfigOption(key, value) 

61 

62 try: 

63 yield 

64 finally: 

65 for key, value in prev.items(): 

66 gdal.SetConfigOption(key, value) 

67 

68 

69def open_raster( 

70 path: str, 

71 mode: str = 'r', 

72 driver: str = '', 

73 default_crs: Optional[gws.Crs] = None, 

74 options: dict = None, 

75) -> 'RasterDataSet': 

76 """Create a raster DataSet from a path. 

77 

78 Args: 

79 path: File path. 

80 mode: 'r' (=read), 'a' (=update), 'w' (=create/write) 

81 driver: Driver name, if omitted, will be suggested from the path extension. 

82 default_crs: Default CRS for geometries (fallback to Webmercator). 

83 options: Options for gdal.OpenEx/CreateDataSource. 

84 """ 

85 

86 dso = _DataSetOptions( 

87 path=path, 

88 mode=mode, 

89 driver=driver, 

90 defaultCrs=default_crs, 

91 gdalOpts=options or {}, 

92 ) 

93 

94 return cast(RasterDataSet, _open(dso, need_raster=True)) 

95 

96 

97def open_vector( 

98 path: str, 

99 mode: str = 'r', 

100 driver: str = '', 

101 encoding: Optional[str] = 'utf8', 

102 default_crs: Optional[gws.Crs] = None, 

103 geometry_as_text: bool = False, 

104 options: dict = None, 

105) -> 'VectorDataSet': 

106 """Create a vector DataSet from a path. 

107 

108 Args: 

109 path: File path. 

110 mode: 'r' (=read), 'a' (=update), 'w' (=create/write) 

111 driver: Driver name, if omitted, will be suggested from the path extension. 

112 encoding: If not None, strings will be automatically decoded. 

113 default_crs: Default CRS for geometries (fallback to Webmercator). 

114 geometry_as_text: Don't interpret geometry, extract raw WKT. 

115 options: Options for gdal.OpenEx/CreateDataSource. 

116 

117 

118 Returns: 

119 DataSet object. 

120 

121 """ 

122 

123 dso = _DataSetOptions( 

124 path=path, 

125 mode=mode, 

126 driver=driver, 

127 defaultCrs=default_crs, 

128 encoding=encoding, 

129 geometryAsText=geometry_as_text, 

130 gdalOpts=options or {}, 

131 ) 

132 

133 return cast(VectorDataSet, _open(dso, need_raster=False)) 

134 

135 

136def open_from_image( 

137 image: gws.Image, 

138 bounds: gws.Bounds, 

139 rotation: gws.Size = None, 

140 options: dict = None, 

141) -> 'RasterDataSet': 

142 """Create an in-memory Dataset from an Image. 

143 

144 Args: 

145 image: Image object. 

146 bounds: Geographic bounds. 

147 x_rotation: GeoTransform x rotation. 

148 y_rotation: GeoTransform y rotation. 

149 options: Driver-specific creation options. 

150 """ 

151 

152 gdal.UseExceptions() 

153 

154 drv = gdal.GetDriverByName('MEM') 

155 img_array = image.to_array() 

156 band_count = img_array.shape[2] 

157 

158 gd = drv.Create( 

159 '', 

160 xsize=img_array.shape[1], 

161 ysize=img_array.shape[0], 

162 bands=band_count, 

163 eType=gdal.GDT_Byte, 

164 options=_option_list(options), 

165 ) 

166 for band in range(band_count): 

167 gd.GetRasterBand(band + 1).WriteArray(img_array[:, :, band]) 

168 

169 gt = _bounds_to_geotransform(bounds, (gd.RasterXSize, gd.RasterYSize), rotation) 

170 

171 gd.SetGeoTransform(gt) 

172 gd.SetSpatialRef(_srs_from_srid(bounds.crs.srid)) 

173 

174 dso = _DataSetOptions(path='') 

175 return RasterDataSet(dso, gd) 

176 

177 

178## 

179 

180 

181class _DriverInfoCache(gws.Data): 

182 infos: list[DriverInfo] 

183 extToName: dict 

184 vectorNames: set[str] 

185 rasterNames: set[str] 

186 

187 

188class _DataSetOptions(gws.Data): 

189 path: str 

190 mode: str 

191 driver: str 

192 encoding: str 

193 defaultCrs: gws.Crs 

194 geometryAsText: bool 

195 gdalOpts: dict 

196 

197 

198def _open(dso: _DataSetOptions, need_raster): 

199 if not dso.mode: 

200 dso.mode = 'r' 

201 if dso.mode not in 'rwa': 

202 raise Error(f'invalid open mode {dso.mode!r}') 

203 

204 gdal.UseExceptions() 

205 

206 drv = _driver_from_args(dso.path, dso.driver, need_raster) 

207 dso.defaultCrs = dso.defaultCrs or gws.lib.crs.WEBMERCATOR 

208 

209 if dso.mode == 'w': 

210 gd = drv.CreateDataSource(dso.path, _option_list(dso.gdalOpts)) 

211 if gd is None: 

212 raise Error(f'cannot create {dso.path!r}') 

213 if need_raster: 

214 return RasterDataSet(dso, gd) 

215 return VectorDataSet(dso, gd) 

216 

217 flags = gdal.OF_VERBOSE_ERROR 

218 if dso.mode == 'r': 

219 flags += gdal.OF_READONLY 

220 if dso.mode == 'a': 

221 flags += gdal.OF_UPDATE 

222 if need_raster: 

223 flags += gdal.OF_RASTER 

224 else: 

225 flags += gdal.OF_VECTOR 

226 

227 gd = gdal.OpenEx(dso.path, flags, open_options=_option_list(dso.gdalOpts)) 

228 if gd is None: 

229 raise Error(f'cannot open {dso.path!r}') 

230 

231 if need_raster: 

232 return RasterDataSet(dso, gd) 

233 return VectorDataSet(dso, gd) 

234 

235 

236class _DataSet: 

237 gdDataset: gdal.Dataset 

238 gdDriver: gdal.Driver 

239 dso: _DataSetOptions 

240 driverName: str 

241 

242 def __init__(self, dso: _DataSetOptions, gd_dataset): 

243 self.gdDataset = gd_dataset 

244 self.gdDriver = self.gdDataset.GetDriver() 

245 self.driverName = self.gdDriver.GetDescription() 

246 self.dso = dso 

247 

248 def __enter__(self): 

249 return self 

250 

251 def __exit__(self, exc_type, exc_val, exc_tb): 

252 self.close() 

253 return False 

254 

255 def close(self): 

256 self.gdDataset.FlushCache() 

257 setattr(self, 'gdDataset', None) 

258 

259 def crs(self) -> Optional[gws.Crs]: 

260 srid = _srid_from_srs(self.gdDataset.GetSpatialRef()) 

261 return gws.lib.crs.get(srid) if srid else None 

262 

263 def set_crs(self, crs: gws.Crs): 

264 srs = _srs_from_srid(crs.srid) 

265 self.gdDataset.SetSpatialRef(srs) 

266 

267 

268class RasterDataSet(_DataSet): 

269 def to_image(self) -> gws.Image: 

270 """Convert the raster dataset to an Image object.""" 

271 

272 band_count = self.gdDataset.RasterCount 

273 x_size = self.gdDataset.RasterXSize 

274 y_size = self.gdDataset.RasterYSize 

275 

276 arr_shape = (y_size, x_size, band_count) 

277 arr = np.zeros(arr_shape, dtype=np.uint8) 

278 

279 for band in range(band_count): 

280 gd_band = self.gdDataset.GetRasterBand(band + 1) 

281 arr[:, :, band] = gd_band.ReadAsArray(0, 0, x_size, y_size) 

282 

283 return gws.lib.image.from_array(arr) 

284 

285 def warp_to_path(self, path: str, options: dict): 

286 """Warp a dataset and store it at the given path. 

287 

288 Args: 

289 path: Destination path. 

290 options: GDAL WarpOptions 

291 

292 See: 

293 https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.WarpOptions 

294 https://gdal.org/en/stable/programs/gdalwarp.html 

295 """ 

296 

297 gdal.UseExceptions() 

298 

299 if 'format' not in options: 

300 options = dict(options) 

301 options['format'] = _driver_from_args(path, '', True).GetDescription() 

302 

303 gd = gdal.Warp(path, self.gdDataset, **options) 

304 if gd is None: 

305 raise Error(f'warp failed') 

306 gd.FlushCache() 

307 gd = None 

308 

309 def save_as(self, path: str, driver: str = '', strict=False, options: dict = None): 

310 """Create a copy of a DataSet. 

311 

312 Args: 

313 path: Destination path. 

314 driver: Driver name, if omitted, will be suggested from the path extension. 

315 strict: If True, fail if some options are not supported. 

316 options: Driver-specific creation options. 

317 """ 

318 

319 gdal.UseExceptions() 

320 

321 drv = _driver_from_args(path, driver, need_raster=True) 

322 gd = drv.CreateCopy( 

323 path, 

324 self.gdDataset, 

325 strict=1 if strict else 0, 

326 options=_option_list(options), 

327 ) 

328 gd.SetMetadata(self.gdDataset.GetMetadata()) 

329 gd.FlushCache() 

330 gd = None 

331 

332 def size(self) -> gws.Size: 

333 return (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize) 

334 

335 def bounds(self) -> gws.Bounds: 

336 return _geotransform_to_bounds( 

337 self.gdDataset.GetGeoTransform(), 

338 (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize), 

339 self.crs() or self.dso.defaultCrs, 

340 ) 

341 

342 

343class VectorDataSet(_DataSet): 

344 @contextlib.contextmanager 

345 def transaction(self): 

346 self.gdDataset.StartTransaction() 

347 try: 

348 yield self 

349 self.gdDataset.CommitTransaction() 

350 except: 

351 self.gdDataset.RollbackTransaction() 

352 raise 

353 

354 def create_layer( 

355 self, 

356 name: str, 

357 columns: dict[str, gws.AttributeType], 

358 geometry_type: gws.GeometryType = None, 

359 crs: gws.Crs = None, 

360 overwrite=False, 

361 options: dict = None, 

362 ) -> 'VectorLayer': 

363 """Create a new layer. 

364 

365 Args: 

366 name: Layer name. 

367 columns: Column definitions. 

368 geometry_type: Geometry type. 

369 crs: CRS for geometries. 

370 overwrite: If True, overwrite existing layer. 

371 options: Driver-specific creation options. 

372 """ 

373 

374 opts = dict(options or {}) 

375 if overwrite: 

376 opts['OVERWRITE'] = 'YES' 

377 enc = (self.dso.encoding or '').upper() 

378 if enc: 

379 driver = self.gdDriver.GetName() 

380 if 'Shapefile' in driver: 

381 opts['ENCODING'] = enc 

382 

383 geom_type = ogr.wkbUnknown 

384 srs = None 

385 

386 if geometry_type: 

387 geom_type = _GEOM_TO_OGR.get(geometry_type) 

388 if not geom_type: 

389 gws.log.warning(f'gdal: unsupported {geometry_type=}') 

390 geom_type = ogr.wkbUnknown 

391 crs = crs or self.dso.defaultCrs 

392 srs = _srs_from_srid(crs.srid) 

393 

394 gd_layer = self.gdDataset.CreateLayer( 

395 name, 

396 geom_type=geom_type, 

397 srs=srs, 

398 options=_option_list(opts), 

399 ) 

400 for col_name, col_type in columns.items(): 

401 fd = ogr.FieldDefn(col_name, _ATTR_TO_OGR[col_type]) 

402 if col_type == gws.AttributeType.bool: 

403 fd.SetSubType(ogr.OFSTBoolean) 

404 gd_layer.CreateField(fd) 

405 

406 return VectorLayer(self, gd_layer) 

407 

408 def layers(self) -> list['VectorLayer']: 

409 """Get all layers.""" 

410 

411 cnt = self.gdDataset.GetLayerCount() 

412 return [VectorLayer(self, self.gdDataset.GetLayerByIndex(n)) for n in range(cnt)] 

413 

414 def layer(self, name_or_index: str | int) -> Optional['VectorLayer']: 

415 """Get a layer by name or index.""" 

416 

417 gd_layer = None 

418 if isinstance(name_or_index, int): 

419 gd_layer = self.gdDataset.GetLayerByIndex(name_or_index) 

420 elif isinstance(name_or_index, str): 

421 gd_layer = self.gdDataset.GetLayerByName(name_or_index) 

422 return VectorLayer(self, gd_layer) if gd_layer else None 

423 

424 def require_layer(self, name_or_index: str | int) -> 'VectorLayer': 

425 """Get a layer by name or index, raise an error if not found.""" 

426 

427 la = self.layer(name_or_index) 

428 if la: 

429 return la 

430 raise Error(f'layer {name_or_index} not found') 

431 

432 

433class VectorLayer: 

434 name: str 

435 dso: _DataSetOptions 

436 gdLayer: ogr.Layer 

437 gdDefn: ogr.FeatureDefn 

438 

439 def __init__(self, ds: VectorDataSet, gd_layer: ogr.Layer): 

440 self.gdLayer = gd_layer 

441 self.gdDefn = self.gdLayer.GetLayerDefn() 

442 self.name = self.gdDefn.GetName() 

443 self.dso = ds.dso 

444 

445 def describe(self) -> gws.DataSetDescription: 

446 desc = gws.DataSetDescription( 

447 columns=[], 

448 columnMap={}, 

449 fullName=self.name, 

450 geometryName='', 

451 geometrySrid=0, 

452 geometryType='', 

453 name=self.name, 

454 schema='', 

455 ) 

456 

457 cols = [] 

458 

459 fid_col = self.gdLayer.GetFIDColumn() 

460 if fid_col: 

461 cols.append( 

462 gws.ColumnDescription( 

463 name=fid_col, 

464 type=_OGR_TO_ATTR[ogr.OFTInteger], 

465 nativeType=ogr.OFTInteger, 

466 isPrimaryKey=True, 

467 columnIndex=0, 

468 ) 

469 ) 

470 

471 for i in range(self.gdDefn.GetFieldCount()): 

472 fdef = self.gdDefn.GetFieldDefn(i) 

473 typ = fdef.GetType() 

474 if typ not in _OGR_TO_ATTR: 

475 continue 

476 attr_type = _OGR_TO_ATTR[typ] 

477 if fdef.GetSubType() == ogr.OFSTBoolean: 

478 attr_type = gws.AttributeType.bool 

479 cols.append( 

480 gws.ColumnDescription( 

481 name=fdef.GetName(), 

482 type=attr_type, 

483 nativeType=typ, 

484 columnIndex=i, 

485 ) 

486 ) 

487 

488 for i in range(self.gdDefn.GetGeomFieldCount()): 

489 fdef = self.gdDefn.GetGeomFieldDefn(i) 

490 typ = fdef.GetType() 

491 cols.append( 

492 gws.ColumnDescription( 

493 name=fdef.GetName() or 'geom', 

494 type=gws.AttributeType.geometry, 

495 nativeType=typ, 

496 columnIndex=i, 

497 geometryType=_OGR_TO_GEOM.get(typ) or gws.GeometryType.geometry, 

498 geometrySrid=_srid_from_srs(fdef.GetSpatialRef()) or self.dso.defaultCrs.srid, 

499 ) 

500 ) 

501 

502 desc.columns = cols 

503 desc.columnMap = {c.name: c for c in cols} 

504 

505 for c in cols: 

506 # NB take the last geom 

507 if c.geometryType: 

508 desc.geometryName = c.name 

509 desc.geometryType = c.geometryType 

510 desc.geometrySrid = c.geometrySrid 

511 

512 return desc 

513 

514 def insert(self, records: list[gws.FeatureRecord]) -> list[int]: 

515 desc = self.describe() 

516 fids = [] 

517 

518 for rec in records: 

519 gd_feature = ogr.Feature(self.gdDefn) 

520 if desc.geometryType and rec.shape: 

521 gd_feature.SetGeometry( 

522 ogr.CreateGeometryFromWkt( 

523 rec.shape.to_wkt(), 

524 _srs_from_srid(rec.shape.crs.srid), 

525 ) 

526 ) 

527 

528 if rec.uid and isinstance(rec.uid, int): 

529 gd_feature.SetFID(rec.uid) 

530 

531 for col in desc.columns: 

532 if col.geometryType or col.isPrimaryKey: 

533 continue 

534 val = rec.attributes.get(col.name) 

535 if val is None: 

536 continue 

537 try: 

538 _attr_to_ogr(gd_feature, int(col.nativeType), col.columnIndex, val, self.dso.encoding) 

539 except Exception as exc: 

540 raise Error(f'field cannot be set: {col.name=} {val=}') from exc 

541 

542 self.gdLayer.CreateFeature(gd_feature) 

543 fids.append(gd_feature.GetFID()) 

544 

545 return fids 

546 

547 def count(self, force=False): 

548 return self.gdLayer.GetFeatureCount(force=1 if force else 0) 

549 

550 def get_all(self) -> list[gws.FeatureRecord]: 

551 return list(self.iter_features()) 

552 

553 def iter_features(self) -> Iterable[gws.FeatureRecord]: 

554 self.gdLayer.ResetReading() 

555 

556 while True: 

557 gd_feature = self.gdLayer.GetNextFeature() 

558 if not gd_feature: 

559 break 

560 yield self._feature_record(gd_feature) 

561 

562 def get(self, fid: int) -> Optional[gws.FeatureRecord]: 

563 gd_feature = self.gdLayer.GetFeature(fid) 

564 if gd_feature: 

565 return self._feature_record(gd_feature) 

566 

567 def _feature_record(self, gd_feature): 

568 rec = gws.FeatureRecord( 

569 attributes={}, 

570 shape=None, 

571 meta={'layerName': self.name}, 

572 uid=str(gd_feature.GetFID()), 

573 ) 

574 

575 for i in range(gd_feature.GetFieldCount()): 

576 fdef = gd_feature.GetFieldDefnRef(i) 

577 val = _attr_from_ogr(gd_feature, fdef.GetType(), fdef.GetSubType(), i, self.dso.encoding) 

578 rec.attributes[fdef.GetName()] = val 

579 

580 cnt = gd_feature.GetGeomFieldCount() 

581 if cnt > 0: 

582 # NB take the last geom 

583 # @TODO multigeometry support 

584 fdef = gd_feature.GetGeomFieldRef(cnt - 1) 

585 if fdef: 

586 srid = _srid_from_srs(fdef.GetSpatialReference()) or self.dso.defaultCrs.srid 

587 if self.dso.geometryAsText: 

588 rec.ewkt = f'SRID={srid};{fdef.ExportToWkt()}' 

589 else: 

590 rec.shape = gws.base.shape.from_wkb(bytes(fdef.ExportToIsoWkb()), gws.lib.crs.get(srid)) 

591 

592 return rec 

593 

594 

595## 

596 

597 

598def _driver_from_args(path, driver_name, need_raster): 

599 di = gws.u.get_app_global('gdal_driver_infos', _fetch_driver_infos) 

600 

601 if not driver_name: 

602 ext = path.split('.')[-1] 

603 names = di.extToName.get(ext) 

604 if not names: 

605 raise Error(f'no default driver found for {path!r}') 

606 if len(names) == 1: 

607 driver_name = names[0] 

608 elif ext in _DEFAULT_DRIVERS: 

609 driver_name = _DEFAULT_DRIVERS[ext] 

610 else: 

611 raise Error(f'multiple drivers found for {path!r}: {names}') 

612 

613 is_vector = driver_name in di.vectorNames 

614 is_raster = driver_name in di.rasterNames 

615 

616 if need_raster: 

617 if not is_raster: 

618 raise Error(f'driver {driver_name!r} is not raster') 

619 return gdal.GetDriverByName(driver_name) 

620 

621 if not is_vector: 

622 raise Error(f'driver {driver_name!r} is not vector') 

623 return ogr.GetDriverByName(driver_name) 

624 

625 

626_DEFAULT_DRIVERS = { 

627 'gif': 'GIF', 

628 'gml': 'GML', 

629 'kml': 'KML', 

630 'tif': 'GTiff', 

631 'tiff': 'GTiff', 

632} 

633 

634 

635def _fetch_driver_infos() -> _DriverInfoCache: 

636 dic = _DriverInfoCache( 

637 infos=[], 

638 extToName={}, 

639 vectorNames=set(), 

640 rasterNames=set(), 

641 ) 

642 

643 for n in range(gdal.GetDriverCount()): 

644 drv = gdal.GetDriver(n) 

645 di = DriverInfo( 

646 index=n, 

647 name=str(drv.ShortName), 

648 longName=str(drv.LongName), 

649 extensions=[], 

650 metaData=dict(drv.GetMetadata() or {}), 

651 ) 

652 dic.infos.append(di) 

653 

654 for e in di.metaData.get(gdal.DMD_EXTENSIONS, '').split(): 

655 dic.extToName.setdefault(e, []).append(di.name) 

656 di.extensions.append(e) 

657 if di.metaData.get('DCAP_VECTOR') == 'YES': 

658 dic.vectorNames.add(di.name) 

659 if di.metaData.get('DCAP_RASTER') == 'YES': 

660 dic.rasterNames.add(di.name) 

661 

662 return dic 

663 

664 

665_name_to_srid = {} 

666 

667 

668def _srs_from_srid(srid): 

669 srs = osr.SpatialReference() 

670 srs.ImportFromEPSG(srid) 

671 return srs 

672 

673 

674def _srid_from_srs(srs): 

675 if not srs: 

676 return 0 

677 

678 name = srs.GetName() 

679 if not name: 

680 wkt = srs.ExportToWkt() 

681 gws.log.warning(f'gdalx: no name for SRS {wkt!r}') 

682 return 0 

683 

684 if name in _name_to_srid: 

685 return _name_to_srid[name] 

686 

687 srid = srs.GetAuthorityCode(None) 

688 if not srid: 

689 wkt = srs.ExportToWkt() 

690 gws.log.warning(f'gdalx: no srid for SRS {wkt!r}') 

691 srid = 0 

692 

693 _name_to_srid[name] = srid 

694 return srid 

695 

696 

697def _attr_from_ogr(gd_feature: ogr.Feature, gtype: int, gsubtype: int, idx: int, encoding: str): 

698 if gd_feature.IsFieldNull(idx): 

699 return None 

700 

701 if gtype == ogr.OFTString: 

702 b = gd_feature.GetFieldAsBinary(idx) 

703 if encoding: 

704 return b.decode(encoding) 

705 return bytes(b) 

706 

707 # GetFieldAsDateTime uses float seconds: 

708 # GetFieldAsDateTime(int i, int *pnYear, int *pnMonth, int *pnDay, int *pnHour, int *pnMinute, float *pfSecond, int *pnTZFlag) 

709 

710 if gtype == ogr.OFTDate: 

711 v = gd_feature.GetFieldAsDateTime(idx) 

712 return datetime.date(v[0], v[1], v[2]) 

713 

714 if gtype == ogr.OFTTime: 

715 v = gd_feature.GetFieldAsDateTime(idx) 

716 sec, fsec = divmod(v[5], 1) 

717 return datetime.time(v[3], v[4], int(sec), int(fsec * 1e6)) 

718 

719 if gtype == ogr.OFTDateTime: 

720 v = gd_feature.GetFieldAsDateTime(idx) 

721 sec, fsec = divmod(v[5], 1) 

722 return datetimex.new(v[0], v[1], v[2], v[3], v[4], int(sec), int(fsec * 1e6), tz=_tzflag_to_tz(v[6])) 

723 

724 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}: 

725 return gd_feature.GetFieldAsIntegerList(idx) 

726 if gtype == ogr.OFTRealList: 

727 return gd_feature.GetFieldAsDoubleList(idx) 

728 if gtype == ogr.OFTStringList: 

729 return list(gd_feature.GetFieldAsStringList(idx)) 

730 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}: 

731 if gsubtype == ogr.OFSTBoolean: 

732 return gd_feature.GetFieldAsInteger(idx) != 0 

733 return gd_feature.GetFieldAsInteger(idx) 

734 if gtype == ogr.OFTReal: 

735 return gd_feature.GetFieldAsDouble(idx) 

736 if gtype == ogr.OFTBinary: 

737 return gd_feature.GetFieldAsBinary(idx) 

738 

739 

740def _tzflag_to_tz(tzflag): 

741 # see gdal/ogr/ogrutils.cpp OGRGetISO8601DateTime 

742 

743 if tzflag == 0 or tzflag == 1: 

744 return '' 

745 if tzflag == 100: 

746 return 'UTC' 

747 if tzflag % 4 != 0: 

748 # @TODO 

749 raise Error(f'unsupported timezone {tzflag=}') 

750 hrs = (100 - tzflag) // 4 

751 return f'Etc/GMT{hrs:+}' 

752 

753 

754def _attr_to_ogr(gd_feature: ogr.Feature, gtype: int, idx: int, value: Any, encoding): 

755 if isinstance(value, decimal.Decimal): 

756 value = float(value) 

757 

758 if gtype == ogr.OFTDate: 

759 return gd_feature.SetField(idx, datetimex.to_iso_date_string(value)) 

760 if gtype == ogr.OFTTime: 

761 return gd_feature.SetField(idx, datetimex.to_iso_time_string(value)) 

762 if gtype == ogr.OFTDateTime: 

763 return gd_feature.SetField(idx, datetimex.to_iso_string(datetimex.to_utc(value), with_tz='Z')) 

764 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}: 

765 return gd_feature.SetField(idx, int(bool(value) if isinstance(value, bool) else value)) 

766 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}: 

767 return gd_feature.SetFieldIntegerList(idx, [int(x) for x in value]) 

768 if gtype == ogr.OFTRealList: 

769 return gd_feature.SetFieldDoubleList(idx, [float(x) for x in value]) 

770 if gtype == ogr.OFTReal: 

771 return gd_feature.SetField(idx, float(value)) 

772 if gtype == ogr.OFTString: 

773 if isinstance(value, bytes): 

774 return gd_feature.SetField(idx, value.decode(encoding or 'utf8')) 

775 return gd_feature.SetField(idx, str(value)) 

776 if gtype == ogr.OFTStringList: 

777 return gd_feature.SetFieldStringList(idx, [str(x) for x in value]) 

778 if gtype == ogr.OFTBinary: 

779 return gd_feature.SetFieldBinaryFromHexString(idx, value.hex() if isinstance(value, bytes) else value) 

780 

781 return gd_feature.SetField(idx, value) 

782 

783 

784def _bounds_to_geotransform(bounds: gws.Bounds, px_size: gws.Size, rotation: gws.Size | None) -> tuple[float, float, float, float, float, float]: 

785 ext = bounds.extent 

786 res_x = (ext[2] - ext[0]) / px_size[0] 

787 res_y = (ext[1] - ext[3]) / px_size[1] 

788 xr = rotation[0] if rotation else 0.0 

789 yr = rotation[1] if rotation else 0.0 

790 return (ext[0], res_x, xr, ext[3], yr, res_y) 

791 

792 

793def _geotransform_to_bounds(gt: tuple[float, float, float, float, float, float], px_size: gws.Size, crs: gws.Crs) -> gws.Bounds: 

794 x0 = gt[0] 

795 x1 = x0 + gt[1] * px_size[0] 

796 y1 = gt[3] 

797 y0 = y1 + gt[5] * px_size[1] 

798 return gws.lib.bounds.from_extent((x0, y0, x1, y1), crs, always_xy=True) 

799 

800 

801def _option_list(opts: dict | None) -> list[str]: 

802 if not opts: 

803 return [] 

804 return [f'{k}={v}' for k, v in opts.items()] 

805 

806 

807_ATTR_TO_OGR = { 

808 gws.AttributeType.bool: ogr.OFTInteger, 

809 gws.AttributeType.bytes: ogr.OFTBinary, 

810 gws.AttributeType.date: ogr.OFTDate, 

811 gws.AttributeType.datetime: ogr.OFTDateTime, 

812 gws.AttributeType.float: ogr.OFTReal, 

813 gws.AttributeType.floatlist: ogr.OFTRealList, 

814 gws.AttributeType.int: ogr.OFTInteger, 

815 gws.AttributeType.intlist: ogr.OFTIntegerList, 

816 gws.AttributeType.str: ogr.OFTString, 

817 gws.AttributeType.strlist: ogr.OFTStringList, 

818 gws.AttributeType.time: ogr.OFTTime, 

819} 

820 

821_OGR_TO_ATTR = { 

822 ogr.OFTBinary: gws.AttributeType.bytes, 

823 ogr.OFTDate: gws.AttributeType.date, 

824 ogr.OFTDateTime: gws.AttributeType.datetime, 

825 ogr.OFTReal: gws.AttributeType.float, 

826 ogr.OFTRealList: gws.AttributeType.floatlist, 

827 ogr.OFTInteger: gws.AttributeType.int, 

828 ogr.OFTIntegerList: gws.AttributeType.intlist, 

829 ogr.OFTInteger64: gws.AttributeType.int, 

830 ogr.OFTInteger64List: gws.AttributeType.intlist, 

831 ogr.OFTString: gws.AttributeType.str, 

832 ogr.OFTStringList: gws.AttributeType.strlist, 

833 ogr.OFTTime: gws.AttributeType.time, 

834} 

835 

836_GEOM_TO_OGR = { 

837 gws.GeometryType.curve: ogr.wkbCurve, 

838 gws.GeometryType.geometrycollection: ogr.wkbGeometryCollection, 

839 gws.GeometryType.linestring: ogr.wkbLineString, 

840 gws.GeometryType.multicurve: ogr.wkbMultiCurve, 

841 gws.GeometryType.multilinestring: ogr.wkbMultiLineString, 

842 gws.GeometryType.multipoint: ogr.wkbMultiPoint, 

843 gws.GeometryType.multipolygon: ogr.wkbMultiPolygon, 

844 gws.GeometryType.multisurface: ogr.wkbMultiSurface, 

845 gws.GeometryType.point: ogr.wkbPoint, 

846 gws.GeometryType.polygon: ogr.wkbPolygon, 

847 gws.GeometryType.polyhedralsurface: ogr.wkbPolyhedralSurface, 

848 gws.GeometryType.surface: ogr.wkbSurface, 

849} 

850 

851_OGR_TO_GEOM = { 

852 ogr.wkbCurve: gws.GeometryType.curve, 

853 ogr.wkbGeometryCollection: gws.GeometryType.geometrycollection, 

854 ogr.wkbLineString: gws.GeometryType.linestring, 

855 ogr.wkbMultiCurve: gws.GeometryType.multicurve, 

856 ogr.wkbMultiLineString: gws.GeometryType.multilinestring, 

857 ogr.wkbMultiPoint: gws.GeometryType.multipoint, 

858 ogr.wkbMultiPolygon: gws.GeometryType.multipolygon, 

859 ogr.wkbMultiSurface: gws.GeometryType.multisurface, 

860 ogr.wkbPoint: gws.GeometryType.point, 

861 ogr.wkbPolygon: gws.GeometryType.polygon, 

862 ogr.wkbPolyhedralSurface: gws.GeometryType.polyhedralsurface, 

863 ogr.wkbSurface: gws.GeometryType.surface, 

864}