Coverage for gws-app/gws/lib/gdalx/__init__.py: 89%
473 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-23 20:06 +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`` or from encoded image bytes with ``open_from_bytes``. 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_driver_infos', _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 driver: str = '',
132 default_crs: Optional[gws.Crs] = None,
133 options: dict = None,
134) -> 'RasterDataSet':
135 """Open a raster data set for reading.
137 Raster data sets are read-only. To write a raster, use ``RasterDataSet.save_as``
138 or ``RasterDataSet.warp_to_path``.
140 Args:
141 path: File path.
142 driver: Driver name. If omitted, the driver is chosen by the path extension.
143 default_crs: CRS to use if the data set has none, web mercator by default.
144 options: Driver-specific open options.
146 Returns:
147 The raster data set.
149 Raises:
150 ``Error``: If no suitable raster driver is found or the data set cannot be opened.
151 """
153 dso = _DataSetOptions(
154 path=path,
155 mode='r',
156 driver=driver,
157 defaultCrs=default_crs,
158 gdalOpts=options or {},
159 )
161 return cast(RasterDataSet, _open(dso, need_raster=True))
164def open_vector(
165 path: str,
166 mode: str = 'r',
167 driver: str = '',
168 encoding: Optional[str] = 'utf8',
169 default_crs: Optional[gws.Crs] = None,
170 geometry_as_text: bool = False,
171 options: dict = None,
172) -> 'VectorDataSet':
173 """Open a vector data set.
175 Args:
176 path: File path.
177 mode: ``r`` (read), ``a`` (update) or ``w`` (create).
178 driver: Driver name. If omitted, the driver is chosen by the path extension.
179 encoding: Encoding of string attributes. If set, strings are decoded when reading,
180 otherwise they are returned as bytes.
181 default_crs: CRS for geometries without one, web mercator by default.
182 geometry_as_text: Do not convert geometries to shapes, return them as EWKT in ``FeatureRecord.ewkt``.
183 options: Driver-specific open or creation options.
185 Returns:
186 The vector data set.
188 Raises:
189 ``Error``: If the mode is invalid, no suitable vector driver is found, or the data set cannot be opened or created.
190 """
192 dso = _DataSetOptions(
193 path=path,
194 mode=mode,
195 driver=driver,
196 defaultCrs=default_crs,
197 encoding=encoding,
198 geometryAsText=geometry_as_text,
199 gdalOpts=options or {},
200 )
202 return cast(VectorDataSet, _open(dso, need_raster=False))
205def open_from_image(
206 image: gws.Image,
207 bounds: gws.Bounds,
208 rotation: gws.Size = None,
209 options: dict = None,
210) -> 'RasterDataSet':
211 """Create an in-memory raster data set from an image.
213 Args:
214 image: Image object.
215 bounds: Bounds of the image.
216 rotation: Geo-transform rotation terms ``(x, y)``, no rotation by default.
217 options: Driver-specific creation options.
219 Returns:
220 The raster data set.
221 """
223 gdal.UseExceptions()
225 drv = gdal.GetDriverByName('MEM')
226 img_array = image.to_array()
227 band_count = img_array.shape[2]
229 gd = drv.Create(
230 '',
231 xsize=img_array.shape[1],
232 ysize=img_array.shape[0],
233 bands=band_count,
234 eType=gdal.GDT_Byte,
235 options=_option_list(options),
236 )
237 for band in range(band_count):
238 gd.GetRasterBand(band + 1).WriteArray(img_array[:, :, band])
240 gt = _bounds_to_geotransform(bounds, (gd.RasterXSize, gd.RasterYSize), rotation)
242 gd.SetGeoTransform(gt)
243 gd.SetSpatialRef(_srs_from_srid(bounds.crs.srid))
245 dso = _DataSetOptions(path='')
246 return RasterDataSet(dso, gd)
249def open_from_bytes(content: bytes, bounds: gws.Bounds) -> 'RasterDataSet':
250 """Create a raster data set from an encoded image.
252 The image is kept in GDAL's memory file system and decoded only when the data set is read.
253 Paletted images are expanded to RGBA.
255 Args:
256 content: Encoded image, in any format GDAL can read.
257 bounds: Bounds of the image.
259 Returns:
260 The raster data set.
262 Raises:
263 ``Error``: If the image cannot be opened.
264 """
266 gdal.UseExceptions()
268 uid = gws.u.random_string(32)
269 src_path = f'/vsimem/{uid}.img'
270 vrt_path = f'/vsimem/{uid}.vrt'
272 gdal.FileFromMemBuffer(src_path, content)
274 try:
275 src = gdal.Open(src_path)
276 except RuntimeError as exc:
277 gdal.Unlink(src_path)
278 raise Error('cannot open image') from exc
280 ext = bounds.extent
281 opts = dict(
282 format='VRT',
283 outputBounds=[ext[0], ext[3], ext[2], ext[1]],
284 outputSRS=bounds.crs.epsg,
285 )
286 if src.GetRasterBand(1).GetColorTable():
287 opts['rgbExpand'] = 'rgba'
289 gd = gdal.Translate(vrt_path, src, **opts)
290 src = None
292 dso = _DataSetOptions(path=vrt_path, memPaths=[src_path, vrt_path])
293 return RasterDataSet(dso, gd)
296##
299class _DriverInfoCache(gws.Data):
300 """Cached information about GDAL drivers."""
302 infos: list[DriverInfo]
303 """All drivers."""
304 extToName: dict
305 """Map of file extensions to lists of driver names."""
306 vectorNames: set[str]
307 """Names of vector drivers."""
308 rasterNames: set[str]
309 """Names of raster drivers."""
312class _DataSetOptions(gws.Data):
313 """Options a data set was opened with."""
315 path: str
316 """File path."""
317 mode: str
318 """Open mode."""
319 driver: str
320 """Driver name."""
321 encoding: str
322 """Encoding of string attributes."""
323 defaultCrs: gws.Crs
324 """CRS to use if the data has none."""
325 geometryAsText: bool
326 """Return geometries as EWKT instead of shapes."""
327 gdalOpts: dict
328 """Driver-specific options."""
329 memPaths: list[str]
330 """Files in GDAL's memory file system, removed on close."""
333def _open(dso: _DataSetOptions, need_raster):
334 if not dso.mode:
335 dso.mode = 'r'
336 if dso.mode not in 'rwa':
337 raise Error(f'invalid open mode {dso.mode!r}')
339 gdal.UseExceptions()
341 drv = _driver_from_args(dso.path, dso.driver, need_raster)
342 dso.defaultCrs = dso.defaultCrs or gws.lib.crs.WEBMERCATOR
344 if dso.mode == 'w':
345 gd = drv.CreateDataSource(dso.path, _option_list(dso.gdalOpts))
346 if gd is None:
347 raise Error(f'cannot create {dso.path!r}')
348 return VectorDataSet(dso, gd)
350 flags = gdal.OF_VERBOSE_ERROR
351 if dso.mode == 'r':
352 flags += gdal.OF_READONLY
353 if dso.mode == 'a':
354 flags += gdal.OF_UPDATE
355 if need_raster:
356 flags += gdal.OF_RASTER
357 else:
358 flags += gdal.OF_VECTOR
360 gd = gdal.OpenEx(dso.path, flags, open_options=_option_list(dso.gdalOpts))
361 if gd is None:
362 raise Error(f'cannot open {dso.path!r}')
364 if need_raster:
365 return RasterDataSet(dso, gd)
366 return VectorDataSet(dso, gd)
369class _DataSet:
370 """Base class for GDAL data sets."""
372 gdDataset: gdal.Dataset
373 """Underlying GDAL data set."""
374 gdDriver: gdal.Driver
375 """Underlying GDAL driver."""
376 dso: _DataSetOptions
377 """Options the data set was opened with."""
378 driverName: str
379 """Driver name."""
381 def __init__(self, dso: _DataSetOptions, gd_dataset):
382 """Wrap a GDAL data set.
384 Args:
385 dso: Options the data set was opened with.
386 gd_dataset: GDAL data set.
387 """
388 self.gdDataset = gd_dataset
389 self.gdDriver = self.gdDataset.GetDriver()
390 self.driverName = self.gdDriver.GetDescription()
391 self.dso = dso
393 def __enter__(self):
394 return self
396 def __exit__(self, exc_type, exc_val, exc_tb):
397 self.close()
398 return False
400 def close(self):
401 """Flush the data set and release it."""
402 self.gdDataset.FlushCache()
403 setattr(self, 'gdDataset', None)
404 for p in self.dso.memPaths or []:
405 gdal.Unlink(p)
407 def crs(self) -> Optional[gws.Crs]:
408 """Get the CRS of the data set.
410 Returns:
411 The CRS, or ``None`` if the data set has no CRS or it is unknown.
412 """
413 srid = _srid_from_srs(self.gdDataset.GetSpatialRef())
414 return gws.lib.crs.get(srid) if srid else None
416 def set_crs(self, crs: gws.Crs):
417 """Set the CRS of the data set.
419 Args:
420 crs: The CRS.
421 """
422 srs = _srs_from_srid(crs.srid)
423 self.gdDataset.SetSpatialRef(srs)
426class RasterDataSet(_DataSet):
427 """Raster data set."""
429 def to_image(self) -> gws.Image:
430 """Convert the raster data set to an image.
432 Each raster band becomes an image channel, values are read as 8-bit.
434 Returns:
435 The image.
436 """
438 band_count = self.gdDataset.RasterCount
439 x_size = self.gdDataset.RasterXSize
440 y_size = self.gdDataset.RasterYSize
442 arr_shape = (y_size, x_size, band_count)
443 arr = np.zeros(arr_shape, dtype=np.uint8)
445 for band in range(band_count):
446 gd_band = self.gdDataset.GetRasterBand(band + 1)
447 arr[:, :, band] = gd_band.ReadAsArray(0, 0, x_size, y_size)
449 return gws.lib.image.from_array(arr)
451 def warp_to_image(self, options: dict) -> gws.Image:
452 """Warp the data set in memory and return the result as an image.
454 See https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.WarpOptions
455 and https://gdal.org/en/stable/programs/gdalwarp.html for the options.
457 Args:
458 options: Keyword arguments for ``gdal.Warp``. The ``format`` option is ignored.
460 Returns:
461 The warped image.
463 Raises:
464 ``Error``: If the warp fails.
465 """
467 gdal.UseExceptions()
469 options = dict(options)
470 options['format'] = 'MEM'
472 gd = gdal.Warp('', self.gdDataset, **options)
473 if gd is None:
474 raise Error(f'warp failed')
476 return RasterDataSet(_DataSetOptions(path=''), gd).to_image()
478 def warp_to_path(self, path: str, options: dict):
479 """Warp the data set and store it at the given path.
481 See https://gdal.org/en/stable/api/python/utilities.html#osgeo.gdal.WarpOptions
482 and https://gdal.org/en/stable/programs/gdalwarp.html for the options.
484 Args:
485 path: Destination path.
486 options: Keyword arguments for ``gdal.Warp``. If ``format`` is not given, it is chosen by the path extension.
488 Raises:
489 ``Error``: If no driver is found for the path or the warp fails.
490 """
492 gdal.UseExceptions()
494 if 'format' not in options:
495 options = dict(options)
496 options['format'] = _driver_from_args(path, '', True).GetDescription()
498 gd = gdal.Warp(path, self.gdDataset, **options)
499 if gd is None:
500 raise Error(f'warp failed')
501 gd.FlushCache()
502 gd = None
504 def save_as(self, path: str, driver: str = '', strict=False, options: dict = None):
505 """Save a copy of the data set, including its metadata.
507 Args:
508 path: Destination path.
509 driver: Driver name. If omitted, the driver is chosen by the path extension.
510 strict: Fail if the copy cannot be made exactly, for example, if the format does not support some data.
511 options: Driver-specific creation options.
513 Raises:
514 ``Error``: If no suitable raster driver is found.
515 """
517 gdal.UseExceptions()
519 drv = _driver_from_args(path, driver, need_raster=True)
520 gd = drv.CreateCopy(
521 path,
522 self.gdDataset,
523 strict=1 if strict else 0,
524 options=_option_list(options),
525 )
526 gd.SetMetadata(self.gdDataset.GetMetadata())
527 gd.FlushCache()
528 gd = None
530 def size(self) -> gws.Size:
531 """Get the raster size.
533 Returns:
534 A ``(width, height)`` tuple in pixels.
535 """
536 return (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize)
538 def bounds(self) -> gws.Bounds:
539 """Get the bounds of the raster, computed from its geo-transform.
541 Returns:
542 The bounds, in the data set CRS or the default CRS.
543 """
544 return _geotransform_to_bounds(
545 self.gdDataset.GetGeoTransform(),
546 (self.gdDataset.RasterXSize, self.gdDataset.RasterYSize),
547 self.crs() or self.dso.defaultCrs,
548 )
551class VectorDataSet(_DataSet):
552 """Vector data set."""
554 @contextlib.contextmanager
555 def transaction(self):
556 """Context manager that runs a transaction.
558 The transaction is committed on success and rolled back on an exception.
560 Yields:
561 This data set.
562 """
563 self.gdDataset.StartTransaction()
564 try:
565 yield self
566 self.gdDataset.CommitTransaction()
567 except:
568 self.gdDataset.RollbackTransaction()
569 raise
571 def create_layer(
572 self,
573 name: str,
574 columns: dict[str, gws.AttributeType],
575 geometry_type: gws.GeometryType = None,
576 crs: gws.Crs = None,
577 overwrite=False,
578 options: dict = None,
579 ) -> 'VectorLayer':
580 """Create a new layer.
582 For Shapefiles, the data set encoding is passed to the driver.
584 Args:
585 name: Layer name.
586 columns: Map of column names to attribute types.
587 geometry_type: Geometry type. If omitted, the layer has no geometry.
588 crs: CRS for geometries, the default CRS of the data set by default.
589 overwrite: Overwrite an existing layer.
590 options: Driver-specific layer creation options.
592 Returns:
593 The new layer.
594 """
596 opts = dict(options or {})
597 if overwrite:
598 opts['OVERWRITE'] = 'YES'
599 enc = (self.dso.encoding or '').upper()
600 if enc:
601 driver = self.gdDriver.GetName()
602 if 'Shapefile' in driver:
603 opts['ENCODING'] = enc
605 geom_type = ogr.wkbUnknown
606 srs = None
608 if geometry_type:
609 geom_type = _GEOM_TO_OGR.get(geometry_type)
610 if not geom_type:
611 gws.log.warning(f'gdal: unsupported {geometry_type=}')
612 geom_type = ogr.wkbUnknown
613 crs = crs or self.dso.defaultCrs
614 srs = _srs_from_srid(crs.srid)
616 gd_layer = self.gdDataset.CreateLayer(
617 name,
618 geom_type=geom_type,
619 srs=srs,
620 options=_option_list(opts),
621 )
622 for col_name, col_type in columns.items():
623 fd = ogr.FieldDefn(col_name, _ATTR_TO_OGR[col_type])
624 if col_type == gws.AttributeType.bool:
625 fd.SetSubType(ogr.OFSTBoolean)
626 gd_layer.CreateField(fd)
628 return VectorLayer(self, gd_layer)
630 def layers(self) -> list['VectorLayer']:
631 """Get all layers.
633 Returns:
634 A list of layers.
635 """
637 cnt = self.gdDataset.GetLayerCount()
638 return [VectorLayer(self, self.gdDataset.GetLayerByIndex(n)) for n in range(cnt)]
640 def layer(self, name_or_index: str | int) -> Optional['VectorLayer']:
641 """Get a layer by name or index.
643 Args:
644 name_or_index: Layer name or index.
646 Returns:
647 The layer, or ``None`` if not found.
648 """
650 gd_layer = None
651 if isinstance(name_or_index, int):
652 gd_layer = self.gdDataset.GetLayerByIndex(name_or_index)
653 elif isinstance(name_or_index, str):
654 gd_layer = self.gdDataset.GetLayerByName(name_or_index)
655 return VectorLayer(self, gd_layer) if gd_layer else None
657 def require_layer(self, name_or_index: str | int) -> 'VectorLayer':
658 """Get a layer by name or index, and fail if it is not found.
660 Args:
661 name_or_index: Layer name or index.
663 Returns:
664 The layer.
666 Raises:
667 ``Error``: If the layer is not found.
668 """
670 la = self.layer(name_or_index)
671 if la:
672 return la
673 raise Error(f'layer {name_or_index} not found')
676class VectorLayer:
677 """Layer of a vector data set."""
679 name: str
680 """Layer name."""
681 dso: _DataSetOptions
682 """Options of the data set."""
683 gdLayer: ogr.Layer
684 """Underlying OGR layer."""
685 gdDefn: ogr.FeatureDefn
686 """Underlying OGR feature definition."""
688 def __init__(self, ds: VectorDataSet, gd_layer: ogr.Layer):
689 """Wrap an OGR layer.
691 Args:
692 ds: Data set the layer belongs to.
693 gd_layer: OGR layer.
694 """
695 self.gdLayer = gd_layer
696 self.gdDefn = self.gdLayer.GetLayerDefn()
697 self.name = self.gdDefn.GetName()
698 self.dso = ds.dso
700 def describe(self) -> gws.DataSetDescription:
701 """Describe the layer columns.
703 The description includes the FID column (as the primary key), attribute columns
704 of supported types, and geometry columns. If there are several geometry columns,
705 the last one is used as the layer geometry.
707 Returns:
708 The layer description.
709 """
710 desc = gws.DataSetDescription(
711 columns=[],
712 columnMap={},
713 fullName=self.name,
714 geometryName='',
715 geometrySrid=0,
716 geometryType='',
717 name=self.name,
718 schema='',
719 )
721 cols = []
723 fid_col = self.gdLayer.GetFIDColumn()
724 if fid_col:
725 cols.append(
726 gws.ColumnDescription(
727 name=fid_col,
728 type=_OGR_TO_ATTR[ogr.OFTInteger],
729 nativeType=ogr.OFTInteger,
730 isPrimaryKey=True,
731 columnIndex=0,
732 )
733 )
735 for i in range(self.gdDefn.GetFieldCount()):
736 fdef = self.gdDefn.GetFieldDefn(i)
737 typ = fdef.GetType()
738 if typ not in _OGR_TO_ATTR:
739 continue
740 attr_type = _OGR_TO_ATTR[typ]
741 if fdef.GetSubType() == ogr.OFSTBoolean:
742 attr_type = gws.AttributeType.bool
743 cols.append(
744 gws.ColumnDescription(
745 name=fdef.GetName(),
746 type=attr_type,
747 nativeType=typ,
748 columnIndex=i,
749 )
750 )
752 for i in range(self.gdDefn.GetGeomFieldCount()):
753 fdef = self.gdDefn.GetGeomFieldDefn(i)
754 typ = fdef.GetType()
755 cols.append(
756 gws.ColumnDescription(
757 name=fdef.GetName() or 'geom',
758 type=gws.AttributeType.geometry,
759 nativeType=typ,
760 columnIndex=i,
761 geometryType=_OGR_TO_GEOM.get(typ) or gws.GeometryType.geometry,
762 geometrySrid=_srid_from_srs(fdef.GetSpatialRef()) or self.dso.defaultCrs.srid,
763 )
764 )
766 desc.columns = cols
767 desc.columnMap = {c.name: c for c in cols}
769 for c in cols:
770 # NB take the last geom
771 if c.geometryType:
772 desc.geometryName = c.name
773 desc.geometryType = c.geometryType
774 desc.geometrySrid = c.geometrySrid
776 return desc
778 def insert(self, records: list[gws.FeatureRecord]) -> list[int]:
779 """Insert feature records into the layer.
781 Integer record uids are used as feature ids. Attributes that are ``None``
782 or have no matching column are skipped.
784 Args:
785 records: Feature records.
787 Returns:
788 Feature ids of the inserted features.
790 Raises:
791 ``Error``: If an attribute value cannot be set.
792 """
793 desc = self.describe()
794 fids = []
796 for rec in records:
797 gd_feature = ogr.Feature(self.gdDefn)
798 if desc.geometryType and rec.shape:
799 gd_feature.SetGeometry(
800 ogr.CreateGeometryFromWkt(
801 rec.shape.to_wkt(),
802 _srs_from_srid(rec.shape.crs.srid),
803 )
804 )
806 if rec.uid and isinstance(rec.uid, int):
807 gd_feature.SetFID(rec.uid)
809 for col in desc.columns:
810 if col.geometryType or col.isPrimaryKey:
811 continue
812 val = rec.attributes.get(col.name)
813 if val is None:
814 continue
815 try:
816 _attr_to_ogr(gd_feature, int(col.nativeType), col.columnIndex, val, self.dso.encoding)
817 except Exception as exc:
818 raise Error(f'field cannot be set: {col.name=} {val=}') from exc
820 self.gdLayer.CreateFeature(gd_feature)
821 fids.append(gd_feature.GetFID())
823 return fids
825 def count(self, force=False):
826 """Count features in the layer.
828 Args:
829 force: Count features even if this is expensive for the driver.
831 Returns:
832 The number of features, or ``-1`` if the count is not available without ``force``.
833 """
834 return self.gdLayer.GetFeatureCount(force=1 if force else 0)
836 def get_all(self) -> list[gws.FeatureRecord]:
837 """Read all features.
839 Returns:
840 A list of feature records.
841 """
842 return list(self.iter_features())
844 def iter_features(self) -> Iterable[gws.FeatureRecord]:
845 """Iterate over all features.
847 Yields:
848 Feature records. The record uid is the feature id as a string,
849 ``meta['layerName']`` is the layer name.
850 """
851 self.gdLayer.ResetReading()
853 while True:
854 gd_feature = self.gdLayer.GetNextFeature()
855 if not gd_feature:
856 break
857 yield self._feature_record(gd_feature)
859 def get(self, fid: int) -> Optional[gws.FeatureRecord]:
860 """Read a feature by its id.
862 Args:
863 fid: Feature id.
865 Returns:
866 The feature record, or ``None`` if not found.
867 """
868 gd_feature = self.gdLayer.GetFeature(fid)
869 if gd_feature:
870 return self._feature_record(gd_feature)
872 def _feature_record(self, gd_feature):
873 rec = gws.FeatureRecord(
874 attributes={},
875 shape=None,
876 meta={'layerName': self.name},
877 uid=str(gd_feature.GetFID()),
878 )
880 for i in range(gd_feature.GetFieldCount()):
881 fdef = gd_feature.GetFieldDefnRef(i)
882 val = _attr_from_ogr(gd_feature, fdef.GetType(), fdef.GetSubType(), i, self.dso.encoding)
883 rec.attributes[fdef.GetName()] = val
885 cnt = gd_feature.GetGeomFieldCount()
886 if cnt > 0:
887 # NB take the last geom
888 # @TODO multigeometry support
889 fdef = gd_feature.GetGeomFieldRef(cnt - 1)
890 if fdef:
891 srid = _srid_from_srs(fdef.GetSpatialReference()) or self.dso.defaultCrs.srid
892 if self.dso.geometryAsText:
893 rec.ewkt = f'SRID={srid};{fdef.ExportToWkt()}'
894 else:
895 rec.shape = gws.lib.shape.from_wkb(bytes(fdef.ExportToIsoWkb()), gws.lib.crs.get(srid))
897 return rec
900##
903def _driver_from_args(path, driver_name, need_raster):
904 di = gws.u.get_app_global('gdal_driver_infos', _fetch_driver_infos)
906 if not driver_name:
907 ext = path.split('.')[-1]
908 names = di.extToName.get(ext)
909 if not names:
910 raise Error(f'no default driver found for {path!r}')
911 if len(names) == 1:
912 driver_name = names[0]
913 elif ext in _DEFAULT_DRIVERS:
914 driver_name = _DEFAULT_DRIVERS[ext]
915 else:
916 raise Error(f'multiple drivers found for {path!r}: {names}')
918 is_vector = driver_name in di.vectorNames
919 is_raster = driver_name in di.rasterNames
921 if need_raster:
922 if not is_raster:
923 raise Error(f'driver {driver_name!r} is not raster')
924 return gdal.GetDriverByName(driver_name)
926 if not is_vector:
927 raise Error(f'driver {driver_name!r} is not vector')
928 return ogr.GetDriverByName(driver_name)
931_DEFAULT_DRIVERS = {
932 'gif': 'GIF',
933 'gml': 'GML',
934 'kml': 'KML',
935 'tif': 'GTiff',
936 'tiff': 'GTiff',
937}
940def _fetch_driver_infos() -> _DriverInfoCache:
941 dic = _DriverInfoCache(
942 infos=[],
943 extToName={},
944 vectorNames=set(),
945 rasterNames=set(),
946 )
948 for n in range(gdal.GetDriverCount()):
949 drv = gdal.GetDriver(n)
950 di = DriverInfo(
951 index=n,
952 name=str(drv.ShortName),
953 longName=str(drv.LongName),
954 extensions=[],
955 metaData=dict(drv.GetMetadata() or {}),
956 )
957 dic.infos.append(di)
959 for e in di.metaData.get(gdal.DMD_EXTENSIONS, '').split():
960 dic.extToName.setdefault(e, []).append(di.name)
961 di.extensions.append(e)
962 if di.metaData.get('DCAP_VECTOR') == 'YES':
963 dic.vectorNames.add(di.name)
964 if di.metaData.get('DCAP_RASTER') == 'YES':
965 dic.rasterNames.add(di.name)
967 return dic
970_name_to_srid = {}
973def _srs_from_srid(srid):
974 # our geotransforms are always x=easting/longitude, so the SRS must use the
975 # traditional axis order, not the authority order (northing/latitude first
976 # for e.g. 4326, 3035, 31466-31469, 3044-3045).
978 srs = osr.SpatialReference()
979 srs.ImportFromEPSG(srid)
980 srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
981 return srs
984def _srid_from_srs(srs):
985 if not srs:
986 return 0
988 name = srs.GetName()
989 if not name:
990 wkt = srs.ExportToWkt()
991 gws.log.warning(f'gdalx: no name for SRS {wkt!r}')
992 return 0
994 if name in _name_to_srid:
995 return _name_to_srid[name]
997 srid = srs.GetAuthorityCode(None)
998 if not srid:
999 wkt = srs.ExportToWkt()
1000 gws.log.warning(f'gdalx: no srid for SRS {wkt!r}')
1001 srid = 0
1003 _name_to_srid[name] = srid
1004 return srid
1007def _attr_from_ogr(gd_feature: ogr.Feature, gtype: int, gsubtype: int, idx: int, encoding: str):
1008 if gd_feature.IsFieldNull(idx):
1009 return None
1011 if gtype == ogr.OFTString:
1012 b = gd_feature.GetFieldAsBinary(idx)
1013 if encoding:
1014 return b.decode(encoding)
1015 return bytes(b)
1017 # GetFieldAsDateTime uses float seconds:
1018 # GetFieldAsDateTime(int i, int *pnYear, int *pnMonth, int *pnDay, int *pnHour, int *pnMinute, float *pfSecond, int *pnTZFlag)
1020 if gtype == ogr.OFTDate:
1021 v = gd_feature.GetFieldAsDateTime(idx)
1022 return datetime.date(v[0], v[1], v[2])
1024 if gtype == ogr.OFTTime:
1025 v = gd_feature.GetFieldAsDateTime(idx)
1026 sec, fsec = divmod(v[5], 1)
1027 return datetime.time(v[3], v[4], int(sec), int(fsec * 1e6))
1029 if gtype == ogr.OFTDateTime:
1030 v = gd_feature.GetFieldAsDateTime(idx)
1031 sec, fsec = divmod(v[5], 1)
1032 return datetimex.new(v[0], v[1], v[2], v[3], v[4], int(sec), int(fsec * 1e6), tz=_tzflag_to_tz(v[6]))
1034 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}:
1035 return gd_feature.GetFieldAsIntegerList(idx)
1036 if gtype == ogr.OFTRealList:
1037 return gd_feature.GetFieldAsDoubleList(idx)
1038 if gtype == ogr.OFTStringList:
1039 return list(gd_feature.GetFieldAsStringList(idx))
1040 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}:
1041 if gsubtype == ogr.OFSTBoolean:
1042 return gd_feature.GetFieldAsInteger(idx) != 0
1043 return gd_feature.GetFieldAsInteger(idx)
1044 if gtype == ogr.OFTReal:
1045 return gd_feature.GetFieldAsDouble(idx)
1046 if gtype == ogr.OFTBinary:
1047 return gd_feature.GetFieldAsBinary(idx)
1050def _tzflag_to_tz(tzflag):
1051 # see gdal/ogr/ogrutils.cpp OGRGetISO8601DateTime
1053 if tzflag == 0 or tzflag == 1:
1054 return ''
1055 if tzflag == 100:
1056 return 'UTC'
1057 if tzflag % 4 != 0:
1058 # @TODO
1059 raise Error(f'unsupported timezone {tzflag=}')
1060 hrs = (100 - tzflag) // 4
1061 return f'Etc/GMT{hrs:+}'
1064def _attr_to_ogr(gd_feature: ogr.Feature, gtype: int, idx: int, value: Any, encoding):
1065 if isinstance(value, decimal.Decimal):
1066 value = float(value)
1068 if gtype == ogr.OFTDate:
1069 return gd_feature.SetField(idx, datetimex.to_iso_date_string(value))
1070 if gtype == ogr.OFTTime:
1071 return gd_feature.SetField(idx, datetimex.to_iso_time_string(value))
1072 if gtype == ogr.OFTDateTime:
1073 return gd_feature.SetField(idx, datetimex.to_iso_string(datetimex.to_utc(value), with_tz='Z'))
1074 if gtype in {ogr.OFTInteger, ogr.OFTInteger64}:
1075 return gd_feature.SetField(idx, int(bool(value) if isinstance(value, bool) else value))
1076 if gtype in {ogr.OFTIntegerList, ogr.OFTInteger64List}:
1077 return gd_feature.SetFieldIntegerList(idx, [int(x) for x in value])
1078 if gtype == ogr.OFTRealList:
1079 return gd_feature.SetFieldDoubleList(idx, [float(x) for x in value])
1080 if gtype == ogr.OFTReal:
1081 return gd_feature.SetField(idx, float(value))
1082 if gtype == ogr.OFTString:
1083 if isinstance(value, bytes):
1084 return gd_feature.SetField(idx, value.decode(encoding or 'utf8'))
1085 return gd_feature.SetField(idx, str(value))
1086 if gtype == ogr.OFTStringList:
1087 return gd_feature.SetFieldStringList(idx, [str(x) for x in value])
1088 if gtype == ogr.OFTBinary:
1089 return gd_feature.SetFieldBinaryFromHexString(idx, value.hex() if isinstance(value, bytes) else value)
1091 return gd_feature.SetField(idx, value)
1094def _bounds_to_geotransform(bounds: gws.Bounds, px_size: gws.Size, rotation: gws.Size | None) -> tuple[float, float, float, float, float, float]:
1095 ext = bounds.extent
1096 res_x = (ext[2] - ext[0]) / px_size[0]
1097 res_y = (ext[1] - ext[3]) / px_size[1]
1098 xr = rotation[0] if rotation else 0.0
1099 yr = rotation[1] if rotation else 0.0
1100 return (ext[0], res_x, xr, ext[3], yr, res_y)
1103def _geotransform_to_bounds(gt: tuple[float, float, float, float, float, float], px_size: gws.Size, crs: gws.Crs) -> gws.Bounds:
1104 x0 = gt[0]
1105 x1 = x0 + gt[1] * px_size[0]
1106 y1 = gt[3]
1107 y0 = y1 + gt[5] * px_size[1]
1108 return gws.lib.bounds.from_extent((x0, y0, x1, y1), crs, always_xy=True)
1111def _option_list(opts: dict | None) -> list[str]:
1112 if not opts:
1113 return []
1114 return [f'{k}={v}' for k, v in opts.items()]
1117_ATTR_TO_OGR = {
1118 gws.AttributeType.bool: ogr.OFTInteger,
1119 gws.AttributeType.bytes: ogr.OFTBinary,
1120 gws.AttributeType.date: ogr.OFTDate,
1121 gws.AttributeType.datetime: ogr.OFTDateTime,
1122 gws.AttributeType.float: ogr.OFTReal,
1123 gws.AttributeType.floatlist: ogr.OFTRealList,
1124 gws.AttributeType.int: ogr.OFTInteger,
1125 gws.AttributeType.intlist: ogr.OFTIntegerList,
1126 gws.AttributeType.str: ogr.OFTString,
1127 gws.AttributeType.strlist: ogr.OFTStringList,
1128 gws.AttributeType.time: ogr.OFTTime,
1129}
1131_OGR_TO_ATTR = {
1132 ogr.OFTBinary: gws.AttributeType.bytes,
1133 ogr.OFTDate: gws.AttributeType.date,
1134 ogr.OFTDateTime: gws.AttributeType.datetime,
1135 ogr.OFTReal: gws.AttributeType.float,
1136 ogr.OFTRealList: gws.AttributeType.floatlist,
1137 ogr.OFTInteger: gws.AttributeType.int,
1138 ogr.OFTIntegerList: gws.AttributeType.intlist,
1139 ogr.OFTInteger64: gws.AttributeType.int,
1140 ogr.OFTInteger64List: gws.AttributeType.intlist,
1141 ogr.OFTString: gws.AttributeType.str,
1142 ogr.OFTStringList: gws.AttributeType.strlist,
1143 ogr.OFTTime: gws.AttributeType.time,
1144}
1146_GEOM_TO_OGR = {
1147 gws.GeometryType.curve: ogr.wkbCurve,
1148 gws.GeometryType.geometrycollection: ogr.wkbGeometryCollection,
1149 gws.GeometryType.linestring: ogr.wkbLineString,
1150 gws.GeometryType.multicurve: ogr.wkbMultiCurve,
1151 gws.GeometryType.multilinestring: ogr.wkbMultiLineString,
1152 gws.GeometryType.multipoint: ogr.wkbMultiPoint,
1153 gws.GeometryType.multipolygon: ogr.wkbMultiPolygon,
1154 gws.GeometryType.multisurface: ogr.wkbMultiSurface,
1155 gws.GeometryType.point: ogr.wkbPoint,
1156 gws.GeometryType.polygon: ogr.wkbPolygon,
1157 gws.GeometryType.polyhedralsurface: ogr.wkbPolyhedralSurface,
1158 gws.GeometryType.surface: ogr.wkbSurface,
1159}
1161_OGR_TO_GEOM = {
1162 ogr.wkbCurve: gws.GeometryType.curve,
1163 ogr.wkbGeometryCollection: gws.GeometryType.geometrycollection,
1164 ogr.wkbLineString: gws.GeometryType.linestring,
1165 ogr.wkbMultiCurve: gws.GeometryType.multicurve,
1166 ogr.wkbMultiLineString: gws.GeometryType.multilinestring,
1167 ogr.wkbMultiPoint: gws.GeometryType.multipoint,
1168 ogr.wkbMultiPolygon: gws.GeometryType.multipolygon,
1169 ogr.wkbMultiSurface: gws.GeometryType.multisurface,
1170 ogr.wkbPoint: gws.GeometryType.point,
1171 ogr.wkbPolygon: gws.GeometryType.polygon,
1172 ogr.wkbPolyhedralSurface: gws.GeometryType.polyhedralsurface,
1173 ogr.wkbSurface: gws.GeometryType.surface,
1174}