Coverage for gws-app/gws/lib/gdalx/__init__.py: 89%
453 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-05 13:35 +0200
1"""GDAL/OGR wrapper.
3This package provides a thin layer over the GDAL Python bindings (``osgeo.gdal``, ``osgeo.ogr``)
4for reading and writing raster and vector data sets.
6Data sets are opened with ``open_raster`` or ``open_vector``, or created in memory
7from an image with ``open_from_image``. If no driver name is given, the driver is chosen
8by the file extension. Data sets are context managers that flush and close themselves on exit.
10- ``RasterDataSet`` reads a raster into a ``gws.Image``, reports its size and bounds,
11 warps it to an image or a file (``gdal.Warp``), and saves copies in other formats.
12- ``VectorDataSet`` gives access to ``VectorLayer`` objects, creates new layers and runs transactions.
13- ``VectorLayer`` describes its columns, reads features as ``gws.FeatureRecord`` objects and inserts records.
15Attribute values are converted between OGR field types and ``gws.AttributeType``,
16geometries between OGR and ``gws.Shape``. Geo-transforms always use the
17easting/longitude-first axis order.
19Example::
21 import gws.lib.gdalx
23 with gws.lib.gdalx.open_raster('/data/ortho.tif') as ds:
24 bounds = ds.bounds()
25 img = ds.to_image()
27 with gws.lib.gdalx.open_vector('/data/out.gpkg', 'w') as ds:
28 la = ds.create_layer('roads', {'name': gws.AttributeType.str}, gws.GeometryType.linestring, crs)
29 la.insert(records)
30"""
32from typing import Any, Optional, Iterable, cast
34import datetime
35import decimal
36import contextlib
37import numpy as np
39from osgeo import gdal
40from osgeo import ogr
41from osgeo import osr
43import gws
44import gws.lib.shape
45import gws.lib.crs
46import gws.lib.bounds
47import gws.lib.image
48import gws.lib.datetimex as datetimex
51class Error(gws.Error):
52 """GDAL error."""
54 pass
57class DriverInfo(gws.Data):
58 """Information about a GDAL driver."""
60 index: int
61 """Driver index in GDAL."""
62 name: str
63 """Short driver name, like ``GTiff``."""
64 longName: str
65 """Long driver name, like ``GeoTIFF``."""
66 extensions: list[str]
67 """File extensions supported by the driver."""
68 metaData: dict
69 """Driver metadata."""
72def get_drivers() -> list[DriverInfo]:
73 """Enumerate GDAL drivers.
75 Returns:
76 Information about all available drivers.
77 """
79 di = gws.u.get_app_global('gdal_drivers', _fetch_driver_infos)
80 return di.infos
83def get_driver(name: str) -> Optional[DriverInfo]:
84 """Get driver info by name.
86 Args:
87 name: Short driver name, like ``GTiff``.
89 Returns:
90 Driver information, or ``None`` if the driver is not found.
91 """
93 for di in get_drivers():
94 if di.name == name:
95 return di
98def supported_attribute_types():
99 """Get attribute types that can be written to vector data sets.
101 Returns:
102 A list of ``gws.AttributeType`` values.
103 """
104 return list(_ATTR_TO_OGR.keys())
107@contextlib.contextmanager
108def gdal_config(options: dict):
109 """Context manager that temporarily sets GDAL config options.
111 The previous values are restored on exit.
113 Args:
114 options: GDAL config options, like ``{'GDAL_CACHEMAX': '512'}``.
115 """
117 prev = {}
118 for key, value in options.items():
119 prev[key] = gdal.GetConfigOption(key)
120 gdal.SetConfigOption(key, value)
122 try:
123 yield
124 finally:
125 for key, value in prev.items():
126 gdal.SetConfigOption(key, value)
129def open_raster(
130 path: str,
131 mode: str = 'r',
132 driver: str = '',
133 default_crs: Optional[gws.Crs] = None,
134 options: dict = None,
135) -> 'RasterDataSet':
136 """Open a raster data set.
138 Args:
139 path: File path.
140 mode: ``r`` (read), ``a`` (update) or ``w`` (create).
141 driver: Driver name. If omitted, the driver is chosen by the path extension.
142 default_crs: CRS to use if the data set has none, web mercator by default.
143 options: Driver-specific open or creation options.
145 Returns:
146 The raster data set.
148 Raises:
149 ``Error``: If the mode is invalid, no suitable raster driver is found, or the data set cannot be opened or created.
150 """
152 dso = _DataSetOptions(
153 path=path,
154 mode=mode,
155 driver=driver,
156 defaultCrs=default_crs,
157 gdalOpts=options or {},
158 )
160 return cast(RasterDataSet, _open(dso, need_raster=True))
163def open_vector(
164 path: str,
165 mode: str = 'r',
166 driver: str = '',
167 encoding: Optional[str] = 'utf8',
168 default_crs: Optional[gws.Crs] = None,
169 geometry_as_text: bool = False,
170 options: dict = None,
171) -> 'VectorDataSet':
172 """Open a vector data set.
174 Args:
175 path: File path.
176 mode: ``r`` (read), ``a`` (update) or ``w`` (create).
177 driver: Driver name. If omitted, the driver is chosen by the path extension.
178 encoding: Encoding of string attributes. If set, strings are decoded when reading,
179 otherwise they are returned as bytes.
180 default_crs: CRS for geometries without one, web mercator by default.
181 geometry_as_text: Do not convert geometries to shapes, return them as EWKT in ``FeatureRecord.ewkt``.
182 options: Driver-specific open or creation options.
184 Returns:
185 The vector data set.
187 Raises:
188 ``Error``: If the mode is invalid, no suitable vector driver is found, or the data set cannot be opened or created.
189 """
191 dso = _DataSetOptions(
192 path=path,
193 mode=mode,
194 driver=driver,
195 defaultCrs=default_crs,
196 encoding=encoding,
197 geometryAsText=geometry_as_text,
198 gdalOpts=options or {},
199 )
201 return cast(VectorDataSet, _open(dso, need_raster=False))
204def open_from_image(
205 image: gws.Image,
206 bounds: gws.Bounds,
207 rotation: gws.Size = None,
208 options: dict = None,
209) -> 'RasterDataSet':
210 """Create an in-memory raster data set from an image.
212 Args:
213 image: Image object.
214 bounds: Bounds of the image.
215 rotation: Geo-transform rotation terms ``(x, y)``, no rotation by default.
216 options: Driver-specific creation options.
218 Returns:
219 The raster data set.
220 """
222 gdal.UseExceptions()
224 drv = gdal.GetDriverByName('MEM')
225 img_array = image.to_array()
226 band_count = img_array.shape[2]
228 gd = drv.Create(
229 '',
230 xsize=img_array.shape[1],
231 ysize=img_array.shape[0],
232 bands=band_count,
233 eType=gdal.GDT_Byte,
234 options=_option_list(options),
235 )
236 for band in range(band_count):
237 gd.GetRasterBand(band + 1).WriteArray(img_array[:, :, band])
239 gt = _bounds_to_geotransform(bounds, (gd.RasterXSize, gd.RasterYSize), rotation)
241 gd.SetGeoTransform(gt)
242 gd.SetSpatialRef(_srs_from_srid(bounds.crs.srid))
244 dso = _DataSetOptions(path='')
245 return RasterDataSet(dso, gd)
248##
251class _DriverInfoCache(gws.Data):
252 """Cached information about GDAL drivers."""
254 infos: list[DriverInfo]
255 """All drivers."""
256 extToName: dict
257 """Map of file extensions to lists of driver names."""
258 vectorNames: set[str]
259 """Names of vector drivers."""
260 rasterNames: set[str]
261 """Names of raster drivers."""
264class _DataSetOptions(gws.Data):
265 """Options a data set was opened with."""
267 path: str
268 """File path."""
269 mode: str
270 """Open mode."""
271 driver: str
272 """Driver name."""
273 encoding: str
274 """Encoding of string attributes."""
275 defaultCrs: gws.Crs
276 """CRS to use if the data has none."""
277 geometryAsText: bool
278 """Return geometries as EWKT instead of shapes."""
279 gdalOpts: dict
280 """Driver-specific options."""
283def _open(dso: _DataSetOptions, need_raster):
284 if not dso.mode:
285 dso.mode = 'r'
286 if dso.mode not in 'rwa':
287 raise Error(f'invalid open mode {dso.mode!r}')
289 gdal.UseExceptions()
291 drv = _driver_from_args(dso.path, dso.driver, need_raster)
292 dso.defaultCrs = dso.defaultCrs or gws.lib.crs.WEBMERCATOR
294 if dso.mode == 'w':
295 gd = drv.CreateDataSource(dso.path, _option_list(dso.gdalOpts))
296 if gd is None:
297 raise Error(f'cannot create {dso.path!r}')
298 if need_raster:
299 return RasterDataSet(dso, gd)
300 return VectorDataSet(dso, gd)
302 flags = gdal.OF_VERBOSE_ERROR
303 if dso.mode == 'r':
304 flags += gdal.OF_READONLY
305 if dso.mode == 'a':
306 flags += gdal.OF_UPDATE
307 if need_raster:
308 flags += gdal.OF_RASTER
309 else:
310 flags += gdal.OF_VECTOR
312 gd = gdal.OpenEx(dso.path, flags, open_options=_option_list(dso.gdalOpts))
313 if gd is None:
314 raise Error(f'cannot open {dso.path!r}')
316 if need_raster:
317 return RasterDataSet(dso, gd)
318 return VectorDataSet(dso, gd)
321class _DataSet:
322 """Base class for GDAL data sets."""
324 gdDataset: gdal.Dataset
325 """Underlying GDAL data set."""
326 gdDriver: gdal.Driver
327 """Underlying GDAL driver."""
328 dso: _DataSetOptions
329 """Options the data set was opened with."""
330 driverName: str
331 """Driver name."""
333 def __init__(self, dso: _DataSetOptions, gd_dataset):
334 """Wrap a GDAL data set.
336 Args:
337 dso: Options the data set was opened with.
338 gd_dataset: GDAL data set.
339 """
340 self.gdDataset = gd_dataset
341 self.gdDriver = self.gdDataset.GetDriver()
342 self.driverName = self.gdDriver.GetDescription()
343 self.dso = dso
345 def __enter__(self):
346 return self
348 def __exit__(self, exc_type, exc_val, exc_tb):
349 self.close()
350 return False
352 def close(self):
353 """Flush the data set and release it."""
354 self.gdDataset.FlushCache()
355 setattr(self, 'gdDataset', None)
357 def crs(self) -> Optional[gws.Crs]:
358 """Get the CRS of the data set.
360 Returns:
361 The CRS, or ``None`` if the data set has no CRS or it is unknown.
362 """
363 srid = _srid_from_srs(self.gdDataset.GetSpatialRef())
364 return gws.lib.crs.get(srid) if srid else None
366 def set_crs(self, crs: gws.Crs):
367 """Set the CRS of the data set.
369 Args:
370 crs: The CRS.
371 """
372 srs = _srs_from_srid(crs.srid)
373 self.gdDataset.SetSpatialRef(srs)
376class RasterDataSet(_DataSet):
377 """Raster data set."""
379 def to_image(self) -> gws.Image:
380 """Convert the raster data set to an image.
382 Each raster band becomes an image channel, values are read as 8-bit.
384 Returns:
385 The image.
386 """
388 band_count = self.gdDataset.RasterCount
389 x_size = self.gdDataset.RasterXSize
390 y_size = self.gdDataset.RasterYSize
392 arr_shape = (y_size, x_size, band_count)
393 arr = np.zeros(arr_shape, dtype=np.uint8)
395 for band in range(band_count):
396 gd_band = self.gdDataset.GetRasterBand(band + 1)
397 arr[:, :, band] = gd_band.ReadAsArray(0, 0, x_size, y_size)
399 return gws.lib.image.from_array(arr)
401 def warp_to_image(self, options: dict) -> gws.Image:
402 """Warp the data set in memory and return the result as an image.
404 See https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.WarpOptions
405 and https://gdal.org/en/stable/programs/gdalwarp.html for the options.
407 Args:
408 options: Keyword arguments for ``gdal.Warp``. The ``format`` option is ignored.
410 Returns:
411 The warped image.
413 Raises:
414 ``Error``: If the warp fails.
415 """
417 gdal.UseExceptions()
419 options = dict(options)
420 options['format'] = 'MEM'
422 gd = gdal.Warp('', self.gdDataset, **options)
423 if gd is None:
424 raise Error(f'warp failed')
426 return RasterDataSet(_DataSetOptions(path=''), gd).to_image()
428 def warp_to_path(self, path: str, options: dict):
429 """Warp the data set and store it at the given path.
431 See https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.WarpOptions
432 and https://gdal.org/en/stable/programs/gdalwarp.html for the options.
434 Args:
435 path: Destination path.
436 options: Keyword arguments for ``gdal.Warp``. If ``format`` is not given, it is chosen by the path extension.
438 Raises:
439 ``Error``: If no driver is found for the path or the warp fails.
440 """
442 gdal.UseExceptions()
444 if 'format' not in options:
445 options = dict(options)
446 options['format'] = _driver_from_args(path, '', True).GetDescription()
448 gd = gdal.Warp(path, self.gdDataset, **options)
449 if gd is None:
450 raise Error(f'warp failed')
451 gd.FlushCache()
452 gd = None
454 def save_as(self, path: str, driver: str = '', strict=False, options: dict = None):
455 """Save a copy of the data set, including its metadata.
457 Args:
458 path: Destination path.
459 driver: Driver name. If omitted, the driver is chosen by the path extension.
460 strict: Fail if the copy cannot be made exactly, for example, if the format does not support some data.
461 options: Driver-specific creation options.
463 Raises:
464 ``Error``: If no suitable raster driver is found.
465 """
467 gdal.UseExceptions()
469 drv = _driver_from_args(path, driver, need_raster=True)
470 gd = drv.CreateCopy(
471 path,
472 self.gdDataset,
473 strict=1 if strict else 0,
474 options=_option_list(options),
475 )
476 gd.SetMetadata(self.gdDataset.GetMetadata())
477 gd.FlushCache()
478 gd = None
480 def size(self) -> gws.Size:
481 """Get the raster size.
483 Returns:
484 A ``(width, height)`` tuple in pixels.
485 """
486 return (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize)
488 def bounds(self) -> gws.Bounds:
489 """Get the bounds of the raster, computed from its geo-transform.
491 Returns:
492 The bounds, in the data set CRS or the default CRS.
493 """
494 return _geotransform_to_bounds(
495 self.gdDataset.GetGeoTransform(),
496 (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize),
497 self.crs() or self.dso.defaultCrs,
498 )
501class VectorDataSet(_DataSet):
502 """Vector data set."""
504 @contextlib.contextmanager
505 def transaction(self):
506 """Context manager that runs a transaction.
508 The transaction is committed on success and rolled back on an exception.
510 Yields:
511 This data set.
512 """
513 self.gdDataset.StartTransaction()
514 try:
515 yield self
516 self.gdDataset.CommitTransaction()
517 except:
518 self.gdDataset.RollbackTransaction()
519 raise
521 def create_layer(
522 self,
523 name: str,
524 columns: dict[str, gws.AttributeType],
525 geometry_type: gws.GeometryType = None,
526 crs: gws.Crs = None,
527 overwrite=False,
528 options: dict = None,
529 ) -> 'VectorLayer':
530 """Create a new layer.
532 For Shapefiles, the data set encoding is passed to the driver.
534 Args:
535 name: Layer name.
536 columns: Map of column names to attribute types.
537 geometry_type: Geometry type. If omitted, the layer has no geometry.
538 crs: CRS for geometries, the default CRS of the data set by default.
539 overwrite: Overwrite an existing layer.
540 options: Driver-specific layer creation options.
542 Returns:
543 The new layer.
544 """
546 opts = dict(options or {})
547 if overwrite:
548 opts['OVERWRITE'] = 'YES'
549 enc = (self.dso.encoding or '').upper()
550 if enc:
551 driver = self.gdDriver.GetName()
552 if 'Shapefile' in driver:
553 opts['ENCODING'] = enc
555 geom_type = ogr.wkbUnknown
556 srs = None
558 if geometry_type:
559 geom_type = _GEOM_TO_OGR.get(geometry_type)
560 if not geom_type:
561 gws.log.warning(f'gdal: unsupported {geometry_type=}')
562 geom_type = ogr.wkbUnknown
563 crs = crs or self.dso.defaultCrs
564 srs = _srs_from_srid(crs.srid)
566 gd_layer = self.gdDataset.CreateLayer(
567 name,
568 geom_type=geom_type,
569 srs=srs,
570 options=_option_list(opts),
571 )
572 for col_name, col_type in columns.items():
573 fd = ogr.FieldDefn(col_name, _ATTR_TO_OGR[col_type])
574 if col_type == gws.AttributeType.bool:
575 fd.SetSubType(ogr.OFSTBoolean)
576 gd_layer.CreateField(fd)
578 return VectorLayer(self, gd_layer)
580 def layers(self) -> list['VectorLayer']:
581 """Get all layers.
583 Returns:
584 A list of layers.
585 """
587 cnt = self.gdDataset.GetLayerCount()
588 return [VectorLayer(self, self.gdDataset.GetLayerByIndex(n)) for n in range(cnt)]
590 def layer(self, name_or_index: str | int) -> Optional['VectorLayer']:
591 """Get a layer by name or index.
593 Args:
594 name_or_index: Layer name or index.
596 Returns:
597 The layer, or ``None`` if not found.
598 """
600 gd_layer = None
601 if isinstance(name_or_index, int):
602 gd_layer = self.gdDataset.GetLayerByIndex(name_or_index)
603 elif isinstance(name_or_index, str):
604 gd_layer = self.gdDataset.GetLayerByName(name_or_index)
605 return VectorLayer(self, gd_layer) if gd_layer else None
607 def require_layer(self, name_or_index: str | int) -> 'VectorLayer':
608 """Get a layer by name or index, and fail if it is not found.
610 Args:
611 name_or_index: Layer name or index.
613 Returns:
614 The layer.
616 Raises:
617 ``Error``: If the layer is not found.
618 """
620 la = self.layer(name_or_index)
621 if la:
622 return la
623 raise Error(f'layer {name_or_index} not found')
626class VectorLayer:
627 """Layer of a vector data set."""
629 name: str
630 """Layer name."""
631 dso: _DataSetOptions
632 """Options of the data set."""
633 gdLayer: ogr.Layer
634 """Underlying OGR layer."""
635 gdDefn: ogr.FeatureDefn
636 """Underlying OGR feature definition."""
638 def __init__(self, ds: VectorDataSet, gd_layer: ogr.Layer):
639 """Wrap an OGR layer.
641 Args:
642 ds: Data set the layer belongs to.
643 gd_layer: OGR layer.
644 """
645 self.gdLayer = gd_layer
646 self.gdDefn = self.gdLayer.GetLayerDefn()
647 self.name = self.gdDefn.GetName()
648 self.dso = ds.dso
650 def describe(self) -> gws.DataSetDescription:
651 """Describe the layer columns.
653 The description includes the FID column (as the primary key), attribute columns
654 of supported types, and geometry columns. If there are several geometry columns,
655 the last one is used as the layer geometry.
657 Returns:
658 The layer description.
659 """
660 desc = gws.DataSetDescription(
661 columns=[],
662 columnMap={},
663 fullName=self.name,
664 geometryName='',
665 geometrySrid=0,
666 geometryType='',
667 name=self.name,
668 schema='',
669 )
671 cols = []
673 fid_col = self.gdLayer.GetFIDColumn()
674 if fid_col:
675 cols.append(
676 gws.ColumnDescription(
677 name=fid_col,
678 type=_OGR_TO_ATTR[ogr.OFTInteger],
679 nativeType=ogr.OFTInteger,
680 isPrimaryKey=True,
681 columnIndex=0,
682 )
683 )
685 for i in range(self.gdDefn.GetFieldCount()):
686 fdef = self.gdDefn.GetFieldDefn(i)
687 typ = fdef.GetType()
688 if typ not in _OGR_TO_ATTR:
689 continue
690 attr_type = _OGR_TO_ATTR[typ]
691 if fdef.GetSubType() == ogr.OFSTBoolean:
692 attr_type = gws.AttributeType.bool
693 cols.append(
694 gws.ColumnDescription(
695 name=fdef.GetName(),
696 type=attr_type,
697 nativeType=typ,
698 columnIndex=i,
699 )
700 )
702 for i in range(self.gdDefn.GetGeomFieldCount()):
703 fdef = self.gdDefn.GetGeomFieldDefn(i)
704 typ = fdef.GetType()
705 cols.append(
706 gws.ColumnDescription(
707 name=fdef.GetName() or 'geom',
708 type=gws.AttributeType.geometry,
709 nativeType=typ,
710 columnIndex=i,
711 geometryType=_OGR_TO_GEOM.get(typ) or gws.GeometryType.geometry,
712 geometrySrid=_srid_from_srs(fdef.GetSpatialRef()) or self.dso.defaultCrs.srid,
713 )
714 )
716 desc.columns = cols
717 desc.columnMap = {c.name: c for c in cols}
719 for c in cols:
720 # NB take the last geom
721 if c.geometryType:
722 desc.geometryName = c.name
723 desc.geometryType = c.geometryType
724 desc.geometrySrid = c.geometrySrid
726 return desc
728 def insert(self, records: list[gws.FeatureRecord]) -> list[int]:
729 """Insert feature records into the layer.
731 Integer record uids are used as feature ids. Attributes that are ``None``
732 or have no matching column are skipped.
734 Args:
735 records: Feature records.
737 Returns:
738 Feature ids of the inserted features.
740 Raises:
741 ``Error``: If an attribute value cannot be set.
742 """
743 desc = self.describe()
744 fids = []
746 for rec in records:
747 gd_feature = ogr.Feature(self.gdDefn)
748 if desc.geometryType and rec.shape:
749 gd_feature.SetGeometry(
750 ogr.CreateGeometryFromWkt(
751 rec.shape.to_wkt(),
752 _srs_from_srid(rec.shape.crs.srid),
753 )
754 )
756 if rec.uid and isinstance(rec.uid, int):
757 gd_feature.SetFID(rec.uid)
759 for col in desc.columns:
760 if col.geometryType or col.isPrimaryKey:
761 continue
762 val = rec.attributes.get(col.name)
763 if val is None:
764 continue
765 try:
766 _attr_to_ogr(gd_feature, int(col.nativeType), col.columnIndex, val, self.dso.encoding)
767 except Exception as exc:
768 raise Error(f'field cannot be set: {col.name=} {val=}') from exc
770 self.gdLayer.CreateFeature(gd_feature)
771 fids.append(gd_feature.GetFID())
773 return fids
775 def count(self, force=False):
776 """Count features in the layer.
778 Args:
779 force: Count features even if this is expensive for the driver.
781 Returns:
782 The number of features, or ``-1`` if the count is not available without ``force``.
783 """
784 return self.gdLayer.GetFeatureCount(force=1 if force else 0)
786 def get_all(self) -> list[gws.FeatureRecord]:
787 """Read all features.
789 Returns:
790 A list of feature records.
791 """
792 return list(self.iter_features())
794 def iter_features(self) -> Iterable[gws.FeatureRecord]:
795 """Iterate over all features.
797 Yields:
798 Feature records. The record uid is the feature id as a string,
799 ``meta['layerName']`` is the layer name.
800 """
801 self.gdLayer.ResetReading()
803 while True:
804 gd_feature = self.gdLayer.GetNextFeature()
805 if not gd_feature:
806 break
807 yield self._feature_record(gd_feature)
809 def get(self, fid: int) -> Optional[gws.FeatureRecord]:
810 """Read a feature by its id.
812 Args:
813 fid: Feature id.
815 Returns:
816 The feature record, or ``None`` if not found.
817 """
818 gd_feature = self.gdLayer.GetFeature(fid)
819 if gd_feature:
820 return self._feature_record(gd_feature)
822 def _feature_record(self, gd_feature):
823 rec = gws.FeatureRecord(
824 attributes={},
825 shape=None,
826 meta={'layerName': self.name},
827 uid=str(gd_feature.GetFID()),
828 )
830 for i in range(gd_feature.GetFieldCount()):
831 fdef = gd_feature.GetFieldDefnRef(i)
832 val = _attr_from_ogr(gd_feature, fdef.GetType(), fdef.GetSubType(), i, self.dso.encoding)
833 rec.attributes[fdef.GetName()] = val
835 cnt = gd_feature.GetGeomFieldCount()
836 if cnt > 0:
837 # NB take the last geom
838 # @TODO multigeometry support
839 fdef = gd_feature.GetGeomFieldRef(cnt - 1)
840 if fdef:
841 srid = _srid_from_srs(fdef.GetSpatialReference()) or self.dso.defaultCrs.srid
842 if self.dso.geometryAsText:
843 rec.ewkt = f'SRID={srid};{fdef.ExportToWkt()}'
844 else:
845 rec.shape = gws.lib.shape.from_wkb(bytes(fdef.ExportToIsoWkb()), gws.lib.crs.get(srid))
847 return rec
850##
853def _driver_from_args(path, driver_name, need_raster):
854 di = gws.u.get_app_global('gdal_driver_infos', _fetch_driver_infos)
856 if not driver_name:
857 ext = path.split('.')[-1]
858 names = di.extToName.get(ext)
859 if not names:
860 raise Error(f'no default driver found for {path!r}')
861 if len(names) == 1:
862 driver_name = names[0]
863 elif ext in _DEFAULT_DRIVERS:
864 driver_name = _DEFAULT_DRIVERS[ext]
865 else:
866 raise Error(f'multiple drivers found for {path!r}: {names}')
868 is_vector = driver_name in di.vectorNames
869 is_raster = driver_name in di.rasterNames
871 if need_raster:
872 if not is_raster:
873 raise Error(f'driver {driver_name!r} is not raster')
874 return gdal.GetDriverByName(driver_name)
876 if not is_vector:
877 raise Error(f'driver {driver_name!r} is not vector')
878 return ogr.GetDriverByName(driver_name)
881_DEFAULT_DRIVERS = {
882 'gif': 'GIF',
883 'gml': 'GML',
884 'kml': 'KML',
885 'tif': 'GTiff',
886 'tiff': 'GTiff',
887}
890def _fetch_driver_infos() -> _DriverInfoCache:
891 dic = _DriverInfoCache(
892 infos=[],
893 extToName={},
894 vectorNames=set(),
895 rasterNames=set(),
896 )
898 for n in range(gdal.GetDriverCount()):
899 drv = gdal.GetDriver(n)
900 di = DriverInfo(
901 index=n,
902 name=str(drv.ShortName),
903 longName=str(drv.LongName),
904 extensions=[],
905 metaData=dict(drv.GetMetadata() or {}),
906 )
907 dic.infos.append(di)
909 for e in di.metaData.get(gdal.DMD_EXTENSIONS, '').split():
910 dic.extToName.setdefault(e, []).append(di.name)
911 di.extensions.append(e)
912 if di.metaData.get('DCAP_VECTOR') == 'YES':
913 dic.vectorNames.add(di.name)
914 if di.metaData.get('DCAP_RASTER') == 'YES':
915 dic.rasterNames.add(di.name)
917 return dic
920_name_to_srid = {}
923def _srs_from_srid(srid):
924 # our geotransforms are always x=easting/longitude, so the SRS must use the
925 # traditional axis order, not the authority order (northing/latitude first
926 # for e.g. 4326, 3035, 31466-31469, 3044-3045).
928 srs = osr.SpatialReference()
929 srs.ImportFromEPSG(srid)
930 srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
931 return srs
934def _srid_from_srs(srs):
935 if not srs:
936 return 0
938 name = srs.GetName()
939 if not name:
940 wkt = srs.ExportToWkt()
941 gws.log.warning(f'gdalx: no name for SRS {wkt!r}')
942 return 0
944 if name in _name_to_srid:
945 return _name_to_srid[name]
947 srid = srs.GetAuthorityCode(None)
948 if not srid:
949 wkt = srs.ExportToWkt()
950 gws.log.warning(f'gdalx: no srid for SRS {wkt!r}')
951 srid = 0
953 _name_to_srid[name] = srid
954 return srid
957def _attr_from_ogr(gd_feature: ogr.Feature, gtype: int, gsubtype: int, idx: int, encoding: str):
958 if gd_feature.IsFieldNull(idx):
959 return None
961 if gtype == ogr.OFTString:
962 b = gd_feature.GetFieldAsBinary(idx)
963 if encoding:
964 return b.decode(encoding)
965 return bytes(b)
967 # GetFieldAsDateTime uses float seconds:
968 # GetFieldAsDateTime(int i, int *pnYear, int *pnMonth, int *pnDay, int *pnHour, int *pnMinute, float *pfSecond, int *pnTZFlag)
970 if gtype == ogr.OFTDate:
971 v = gd_feature.GetFieldAsDateTime(idx)
972 return datetime.date(v[0], v[1], v[2])
974 if gtype == ogr.OFTTime:
975 v = gd_feature.GetFieldAsDateTime(idx)
976 sec, fsec = divmod(v[5], 1)
977 return datetime.time(v[3], v[4], int(sec), int(fsec * 1e6))
979 if gtype == ogr.OFTDateTime:
980 v = gd_feature.GetFieldAsDateTime(idx)
981 sec, fsec = divmod(v[5], 1)
982 return datetimex.new(v[0], v[1], v[2], v[3], v[4], int(sec), int(fsec * 1e6), tz=_tzflag_to_tz(v[6]))
984 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}:
985 return gd_feature.GetFieldAsIntegerList(idx)
986 if gtype == ogr.OFTRealList:
987 return gd_feature.GetFieldAsDoubleList(idx)
988 if gtype == ogr.OFTStringList:
989 return list(gd_feature.GetFieldAsStringList(idx))
990 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}:
991 if gsubtype == ogr.OFSTBoolean:
992 return gd_feature.GetFieldAsInteger(idx) != 0
993 return gd_feature.GetFieldAsInteger(idx)
994 if gtype == ogr.OFTReal:
995 return gd_feature.GetFieldAsDouble(idx)
996 if gtype == ogr.OFTBinary:
997 return gd_feature.GetFieldAsBinary(idx)
1000def _tzflag_to_tz(tzflag):
1001 # see gdal/ogr/ogrutils.cpp OGRGetISO8601DateTime
1003 if tzflag == 0 or tzflag == 1:
1004 return ''
1005 if tzflag == 100:
1006 return 'UTC'
1007 if tzflag % 4 != 0:
1008 # @TODO
1009 raise Error(f'unsupported timezone {tzflag=}')
1010 hrs = (100 - tzflag) // 4
1011 return f'Etc/GMT{hrs:+}'
1014def _attr_to_ogr(gd_feature: ogr.Feature, gtype: int, idx: int, value: Any, encoding):
1015 if isinstance(value, decimal.Decimal):
1016 value = float(value)
1018 if gtype == ogr.OFTDate:
1019 return gd_feature.SetField(idx, datetimex.to_iso_date_string(value))
1020 if gtype == ogr.OFTTime:
1021 return gd_feature.SetField(idx, datetimex.to_iso_time_string(value))
1022 if gtype == ogr.OFTDateTime:
1023 return gd_feature.SetField(idx, datetimex.to_iso_string(datetimex.to_utc(value), with_tz='Z'))
1024 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}:
1025 return gd_feature.SetField(idx, int(bool(value) if isinstance(value, bool) else value))
1026 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}:
1027 return gd_feature.SetFieldIntegerList(idx, [int(x) for x in value])
1028 if gtype == ogr.OFTRealList:
1029 return gd_feature.SetFieldDoubleList(idx, [float(x) for x in value])
1030 if gtype == ogr.OFTReal:
1031 return gd_feature.SetField(idx, float(value))
1032 if gtype == ogr.OFTString:
1033 if isinstance(value, bytes):
1034 return gd_feature.SetField(idx, value.decode(encoding or 'utf8'))
1035 return gd_feature.SetField(idx, str(value))
1036 if gtype == ogr.OFTStringList:
1037 return gd_feature.SetFieldStringList(idx, [str(x) for x in value])
1038 if gtype == ogr.OFTBinary:
1039 return gd_feature.SetFieldBinaryFromHexString(idx, value.hex() if isinstance(value, bytes) else value)
1041 return gd_feature.SetField(idx, value)
1044def _bounds_to_geotransform(bounds: gws.Bounds, px_size: gws.Size, rotation: gws.Size | None) -> tuple[float, float, float, float, float, float]:
1045 ext = bounds.extent
1046 res_x = (ext[2] - ext[0]) / px_size[0]
1047 res_y = (ext[1] - ext[3]) / px_size[1]
1048 xr = rotation[0] if rotation else 0.0
1049 yr = rotation[1] if rotation else 0.0
1050 return (ext[0], res_x, xr, ext[3], yr, res_y)
1053def _geotransform_to_bounds(gt: tuple[float, float, float, float, float, float], px_size: gws.Size, crs: gws.Crs) -> gws.Bounds:
1054 x0 = gt[0]
1055 x1 = x0 + gt[1] * px_size[0]
1056 y1 = gt[3]
1057 y0 = y1 + gt[5] * px_size[1]
1058 return gws.lib.bounds.from_extent((x0, y0, x1, y1), crs, always_xy=True)
1061def _option_list(opts: dict | None) -> list[str]:
1062 if not opts:
1063 return []
1064 return [f'{k}={v}' for k, v in opts.items()]
1067_ATTR_TO_OGR = {
1068 gws.AttributeType.bool: ogr.OFTInteger,
1069 gws.AttributeType.bytes: ogr.OFTBinary,
1070 gws.AttributeType.date: ogr.OFTDate,
1071 gws.AttributeType.datetime: ogr.OFTDateTime,
1072 gws.AttributeType.float: ogr.OFTReal,
1073 gws.AttributeType.floatlist: ogr.OFTRealList,
1074 gws.AttributeType.int: ogr.OFTInteger,
1075 gws.AttributeType.intlist: ogr.OFTIntegerList,
1076 gws.AttributeType.str: ogr.OFTString,
1077 gws.AttributeType.strlist: ogr.OFTStringList,
1078 gws.AttributeType.time: ogr.OFTTime,
1079}
1081_OGR_TO_ATTR = {
1082 ogr.OFTBinary: gws.AttributeType.bytes,
1083 ogr.OFTDate: gws.AttributeType.date,
1084 ogr.OFTDateTime: gws.AttributeType.datetime,
1085 ogr.OFTReal: gws.AttributeType.float,
1086 ogr.OFTRealList: gws.AttributeType.floatlist,
1087 ogr.OFTInteger: gws.AttributeType.int,
1088 ogr.OFTIntegerList: gws.AttributeType.intlist,
1089 ogr.OFTInteger64: gws.AttributeType.int,
1090 ogr.OFTInteger64List: gws.AttributeType.intlist,
1091 ogr.OFTString: gws.AttributeType.str,
1092 ogr.OFTStringList: gws.AttributeType.strlist,
1093 ogr.OFTTime: gws.AttributeType.time,
1094}
1096_GEOM_TO_OGR = {
1097 gws.GeometryType.curve: ogr.wkbCurve,
1098 gws.GeometryType.geometrycollection: ogr.wkbGeometryCollection,
1099 gws.GeometryType.linestring: ogr.wkbLineString,
1100 gws.GeometryType.multicurve: ogr.wkbMultiCurve,
1101 gws.GeometryType.multilinestring: ogr.wkbMultiLineString,
1102 gws.GeometryType.multipoint: ogr.wkbMultiPoint,
1103 gws.GeometryType.multipolygon: ogr.wkbMultiPolygon,
1104 gws.GeometryType.multisurface: ogr.wkbMultiSurface,
1105 gws.GeometryType.point: ogr.wkbPoint,
1106 gws.GeometryType.polygon: ogr.wkbPolygon,
1107 gws.GeometryType.polyhedralsurface: ogr.wkbPolyhedralSurface,
1108 gws.GeometryType.surface: ogr.wkbSurface,
1109}
1111_OGR_TO_GEOM = {
1112 ogr.wkbCurve: gws.GeometryType.curve,
1113 ogr.wkbGeometryCollection: gws.GeometryType.geometrycollection,
1114 ogr.wkbLineString: gws.GeometryType.linestring,
1115 ogr.wkbMultiCurve: gws.GeometryType.multicurve,
1116 ogr.wkbMultiLineString: gws.GeometryType.multilinestring,
1117 ogr.wkbMultiPoint: gws.GeometryType.multipoint,
1118 ogr.wkbMultiPolygon: gws.GeometryType.multipolygon,
1119 ogr.wkbMultiSurface: gws.GeometryType.multisurface,
1120 ogr.wkbPoint: gws.GeometryType.point,
1121 ogr.wkbPolygon: gws.GeometryType.polygon,
1122 ogr.wkbPolyhedralSurface: gws.GeometryType.polyhedralsurface,
1123 ogr.wkbSurface: gws.GeometryType.surface,
1124}