Tip

For an interactive online version click here:

Example: Working with geodatasets

from hydromt import DataCatalog
from hydromt.log import setuplog

logger = setuplog("raster data", log_level=10)
dc = DataCatalog(logger=logger, data_libs=["artifact_data"])

2024-04-25 17:20:29,751 - raster data - log - INFO - HydroMT version: 0.9.5.dev0

2024-04-25 17:20:29,776 - raster data - data_catalog - INFO - Reading data catalog archive artifact_data v0.0.8

2024-04-25 17:20:29,777 - raster data - data_catalog - INFO - Parsing data catalog from /home/runner/.hydromt_data/artifact_data/v0.0.8/data_catalog.yml

Here, we illustrate some common GIS problems and how the functionality of the DataArray/Dataset vector accessor can be used. The data is accessed using the HydroMT DataCatalog. For more information see the Reading vector data example.

Geospatial attributes

Some of the available geospatial attributes are listed here below. For all of attributes (and methods for that matter): check the HydroMT API reference

# Get the waterlevel dataset from the datacatalog
# The waterlevels are series per location (points; here in lat lon)
ds = dc.get_geodataset("gtsmv3_eu_era5")

# Coordinate reference system
ds.vector.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# Geospatial data as geometry objects (we show the first 5 points)
ds.vector.geometry.head()

stations
13670    POINT (12.25342 45.25635)
2798     POINT (12.22412 45.27100)
2799     POINT (12.29736 45.22705)
13775    POINT (12.75879 45.24902)
13723    POINT (12.50244 45.25635)
dtype: geometry

# names of x- and y coordinates

(ds.vector.x_name, ds.vector.y_name)

('lon', 'lat')

# Bounding box of the geospatial data
ds.vector.bounds

array([12.22412, 45.22705, 12.99316, 45.62256])

Reprojection

# Reproject data to Pseudo Mercator (EPSG: 3857)
ds_pm = ds.vector.to_crs(3857)

# Coordinate reference system
ds_pm.vector.crs

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# Check that the coordinate values are indeed no longer in degrees
ds_pm.lon.values[0]

1364044.4748761142

Conversion

# Create an ogr compliant dataset from ds
# When written in netcdf4 format, this can be read by ogr (osgeo; QGIS)
from numpy import mean

ds_ogr = ds.vector.ogr_compliant(reducer=mean)
ds_ogr

<xarray.Dataset> Size: 388B
Dimensions:      (stations: 19)
Coordinates:
  * stations     (stations) int32 76B 13670 2798 2799 13775 ... 2791 2790 2789
    ogc_wkt      (stations) object 152B 'POINT (12.25342 45.25635)' ... 'POIN...
    spatial_ref  int64 8B 0
Data variables:
    waterlevel   (stations) float64 152B dask.array<chunksize=(19,), meta=np.ndarray>
Attributes:
    Conventions:         CF-1.6
    GDAL:                GDAL 3.8.4
    ogr_geometry_field:  ogc_wkt
    ogr_layer_type:      Point

xarray.Dataset

• Dimensions:
  • stations: 19
• Coordinates: (3)
  • stations
    (stations)
    int32
    13670 2798 2799 ... 2791 2790 2789

    array([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype=int32)

  • ogc_wkt
    (stations)
    object
    'POINT (12.25342 45.25635)' ... ...

    long_name :

    Geometry as ISO WKT

    grid_mapping :

    spatial_ref

    array(['POINT (12.25342 45.25635)', 'POINT (12.22412 45.271)',
           'POINT (12.29736 45.22705)', 'POINT (12.75879 45.24902)',
           'POINT (12.50244 45.25635)', 'POINT (12.26807 45.34424)',
           'POINT (12.99316 45.24902)', 'POINT (12.34131 45.37353)',
           'POINT (12.4585 45.41748)', 'POINT (12.4292 45.41748)',
           'POINT (12.4292 45.41748)', 'POINT (12.50244 45.46143)',
           'POINT (12.59033 45.49072)', 'POINT (12.75146 45.50537)',
           'POINT (12.50244 45.50537)', 'POINT (12.67822 45.52002)',
           'POINT (12.76611 45.54932)', 'POINT (12.83936 45.57861)',
           'POINT (12.92725 45.62256)'], dtype=object)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    array(0)

• Data variables: (1)
  • waterlevel
    (stations)
    float64
    dask.array<chunksize=(19,), meta=np.ndarray>

    ogr_field_name :

    waterlevel

    ogr_field_type :

    Real

    Array Chunk Bytes 152 B 152 B Shape (19,) (19,) Dask graph 1 chunks in 6 graph layers Data type float64 numpy.ndarray

• Indexes: (1)
  • stations
    PandasIndex

    PandasIndex(Index([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype='int32', name='stations'))

• Attributes: (4)

  Conventions :

  CF-1.6

  GDAL :

  GDAL 3.8.4

  ogr_geometry_field :

  ogc_wkt

  ogr_layer_type :

  Point

# Convert the geospatial data to wkt strings
ds_wkt = ds.vector.to_wkt()
ds_wkt

<xarray.DataArray 'waterlevel' (time: 2016, stations: 19)> Size: 306kB
dask.array<open_dataset-waterlevel, shape=(2016, 19), dtype=float64, chunksize=(2016, 19), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 16kB 2010-02-01 ... 2010-02-14T23:50:00
    spatial_ref  int32 4B ...
  * stations     (stations) int32 76B 13670 2798 2799 13775 ... 2791 2790 2789
    ogc_wkt      (stations) object 152B 'POINT (12.25342 45.25635)' ... 'POIN...
Attributes:
    long_name:       sea_surface_height_above_mean_sea_level
    units:           m
    CDI_grid_type:   unstructured
    short_name:      waterlevel
    description:     Total water level resulting from the combination of baro...
    category:        ocean
    paper_doi:       10.24381/cds.8c59054f
    paper_ref:       Copernicus Climate Change Service 2019
    source_license:  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licenc...
    source_url:      https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24...
    source_version:  GTSM v3.0

xarray.DataArray

'waterlevel'

• time: 2016
• stations: 19

• dask.array<chunksize=(2016, 19), meta=np.ndarray>

  Array Chunk Bytes 299.25 kiB 299.25 kiB Shape (2016, 19) (2016, 19) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2010-02-01 ... 2010-02-14T23:50:00

    axis :

    T

    long_name :

    time

    short_name :

    time

    array(['2010-02-01T00:00:00.000000000', '2010-02-01T00:10:00.000000000',
           '2010-02-01T00:20:00.000000000', ..., '2010-02-14T23:30:00.000000000',
           '2010-02-14T23:40:00.000000000', '2010-02-14T23:50:00.000000000'],
          dtype='datetime64[ns]')

  • spatial_ref
    ()
    int32
    ...

    x_dim :

    lon

    y_dim :

    lat

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    geom_format :

    wkt

    index_dim :

    stations

    geom_name :

    ogc_wkt

    [1 values with dtype=int32]

  • stations
    (stations)
    int32
    13670 2798 2799 ... 2791 2790 2789

    array([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype=int32)

  • ogc_wkt
    (stations)
    object
    'POINT (12.25342 45.25635)' ... ...

    array(['POINT (12.25342 45.25635)', 'POINT (12.22412 45.271)',
           'POINT (12.29736 45.22705)', 'POINT (12.75879 45.24902)',
           'POINT (12.50244 45.25635)', 'POINT (12.26807 45.34424)',
           'POINT (12.99316 45.24902)', 'POINT (12.34131 45.37353)',
           'POINT (12.4585 45.41748)', 'POINT (12.4292 45.41748)',
           'POINT (12.4292 45.41748)', 'POINT (12.50244 45.46143)',
           'POINT (12.59033 45.49072)', 'POINT (12.75146 45.50537)',
           'POINT (12.50244 45.50537)', 'POINT (12.67822 45.52002)',
           'POINT (12.76611 45.54932)', 'POINT (12.83936 45.57861)',
           'POINT (12.92725 45.62256)'], dtype=object)

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-02-01 00:00:00', '2010-02-01 00:10:00',
                   '2010-02-01 00:20:00', '2010-02-01 00:30:00',
                   '2010-02-01 00:40:00', '2010-02-01 00:50:00',
                   '2010-02-01 01:00:00', '2010-02-01 01:10:00',
                   '2010-02-01 01:20:00', '2010-02-01 01:30:00',
                   ...
                   '2010-02-14 22:20:00', '2010-02-14 22:30:00',
                   '2010-02-14 22:40:00', '2010-02-14 22:50:00',
                   '2010-02-14 23:00:00', '2010-02-14 23:10:00',
                   '2010-02-14 23:20:00', '2010-02-14 23:30:00',
                   '2010-02-14 23:40:00', '2010-02-14 23:50:00'],
                  dtype='datetime64[ns]', name='time', length=2016, freq=None))

  • stations
    PandasIndex

    PandasIndex(Index([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype='int32', name='stations'))

• Attributes: (11)

  long_name :

  sea_surface_height_above_mean_sea_level

  units :

  m

  CDI_grid_type :

  unstructured

  short_name :

  waterlevel

  description :

  Total water level resulting from the combination of barotropic tides and surges and mean sea-level

  category :

  ocean

  paper_doi :

  10.24381/cds.8c59054f

  paper_ref :

  Copernicus Climate Change Service 2019

  source_license :

  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licence-to-use-copernicus-products

  source_url :

  https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24381/cds.8c59054f?tab=overview

  source_version :

  GTSM v3.0

# Convert the geospatial data to geometry objects
ds_geom = ds.vector.to_geom()
ds_geom

<xarray.DataArray 'waterlevel' (time: 2016, stations: 19)> Size: 306kB
dask.array<open_dataset-waterlevel, shape=(2016, 19), dtype=float64, chunksize=(2016, 19), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 16kB 2010-02-01 ... 2010-02-14T23:50:00
    spatial_ref  int32 4B ...
  * stations     (stations) int32 76B 13670 2798 2799 13775 ... 2791 2790 2789
    geometry     (stations) object 152B POINT (12.25342 45.25635) ... POINT (...
Attributes:
    long_name:       sea_surface_height_above_mean_sea_level
    units:           m
    CDI_grid_type:   unstructured
    short_name:      waterlevel
    description:     Total water level resulting from the combination of baro...
    category:        ocean
    paper_doi:       10.24381/cds.8c59054f
    paper_ref:       Copernicus Climate Change Service 2019
    source_license:  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licenc...
    source_url:      https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24...
    source_version:  GTSM v3.0

xarray.DataArray

'waterlevel'

• time: 2016
• stations: 19

• dask.array<chunksize=(2016, 19), meta=np.ndarray>

  Array Chunk Bytes 299.25 kiB 299.25 kiB Shape (2016, 19) (2016, 19) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2010-02-01 ... 2010-02-14T23:50:00

    axis :

    T

    long_name :

    time

    short_name :

    time

    array(['2010-02-01T00:00:00.000000000', '2010-02-01T00:10:00.000000000',
           '2010-02-01T00:20:00.000000000', ..., '2010-02-14T23:30:00.000000000',
           '2010-02-14T23:40:00.000000000', '2010-02-14T23:50:00.000000000'],
          dtype='datetime64[ns]')

  • spatial_ref
    ()
    int32
    ...

    x_dim :

    lon

    y_dim :

    lat

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    geom_format :

    geom

    index_dim :

    stations

    geom_name :

    geometry

    [1 values with dtype=int32]

  • stations
    (stations)
    int32
    13670 2798 2799 ... 2791 2790 2789

    array([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype=int32)

  • geometry
    (stations)
    object
    POINT (12.25342 45.25635) ... PO...

    array([<POINT (12.253 45.256)>, <POINT (12.224 45.271)>,
           <POINT (12.297 45.227)>, <POINT (12.759 45.249)>,
           <POINT (12.502 45.256)>, <POINT (12.268 45.344)>,
           <POINT (12.993 45.249)>, <POINT (12.341 45.374)>,
           <POINT (12.458 45.417)>, <POINT (12.429 45.417)>,
           <POINT (12.429 45.417)>, <POINT (12.502 45.461)>,
           <POINT (12.59 45.491)>, <POINT (12.751 45.505)>,
           <POINT (12.502 45.505)>, <POINT (12.678 45.52)>,
           <POINT (12.766 45.549)>, <POINT (12.839 45.579)>,
           <POINT (12.927 45.623)>], dtype=object)

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-02-01 00:00:00', '2010-02-01 00:10:00',
                   '2010-02-01 00:20:00', '2010-02-01 00:30:00',
                   '2010-02-01 00:40:00', '2010-02-01 00:50:00',
                   '2010-02-01 01:00:00', '2010-02-01 01:10:00',
                   '2010-02-01 01:20:00', '2010-02-01 01:30:00',
                   ...
                   '2010-02-14 22:20:00', '2010-02-14 22:30:00',
                   '2010-02-14 22:40:00', '2010-02-14 22:50:00',
                   '2010-02-14 23:00:00', '2010-02-14 23:10:00',
                   '2010-02-14 23:20:00', '2010-02-14 23:30:00',
                   '2010-02-14 23:40:00', '2010-02-14 23:50:00'],
                  dtype='datetime64[ns]', name='time', length=2016, freq=None))

  • stations
    PandasIndex

    PandasIndex(Index([13670,  2798,  2799, 13775, 13723,  2797, 13822,  2796, 23305, 22721,
            2795,  2794,  2793, 13774, 13722,  2792,  2791,  2790,  2789],
          dtype='int32', name='stations'))

• Attributes: (11)

  long_name :

  sea_surface_height_above_mean_sea_level

  units :

  m

  CDI_grid_type :

  unstructured

  short_name :

  waterlevel

  description :

  Total water level resulting from the combination of barotropic tides and surges and mean sea-level

  category :

  ocean

  paper_doi :

  10.24381/cds.8c59054f

  paper_ref :

  Copernicus Climate Change Service 2019

  source_license :

  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licence-to-use-copernicus-products

  source_url :

  https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24381/cds.8c59054f?tab=overview

  source_version :

  GTSM v3.0

I/O (Internal and External)

# Create a dummy GeoDataFrame for this example
from geopandas import GeoDataFrame
from shapely.geometry import Point
from pyproj import CRS
from hydromt.vector import GeoDataset, GeoDataArray
import numpy as np

gdf = GeoDataFrame(
    [{"Loc": "I", "Stuff": 1}, {"Loc": "II", "Stuff": 2}],
    geometry=[Point(0, 0), Point(1, 1)],
    crs=CRS.from_epsg(4326),
)

ds_gdf = GeoDataArray.from_gdf(gdf, np.arange(gdf.index.size))
ds_gdf

<xarray.DataArray (index: 2)> Size: 16B
array([0, 1])
Coordinates:
  * index        (index) int64 16B 0 1
    Loc          (index) object 16B 'I' 'II'
    Stuff        (index) int64 16B 1 2
    geometry     (index) object 16B POINT (0 0) POINT (1 1)
    spatial_ref  int64 8B 0

xarray.DataArray

• index: 2

• 0 1

  array([0, 1])

• Coordinates: (5)
  • index
    (index)
    int64
    0 1

    array([0, 1])

  • Loc
    (index)
    object
    'I' 'II'

    array(['I', 'II'], dtype=object)

  • Stuff
    (index)
    int64
    1 2

    array([1, 2])

  • geometry
    (index)
    object
    POINT (0 0) POINT (1 1)

    array([<POINT (0 0)>, <POINT (1 1)>], dtype=object)

  • spatial_ref
    ()
    int64
    0

    geom_name :

    geometry

    geom_format :

    geom

    index_dim :

    index

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    array(0)

• Indexes: (1)
  • index
    PandasIndex

    PandasIndex(Index([0, 1], dtype='int64', name='index'))

• Attributes: (0)

# Write the Dataset of the DataCatalog to a GeoDataFrame
# Waterlevel has besides stations a time dimension
# GeoDataFrames don't like vectors/ 2d array's, so if you want to keep the variable it can be reduced along the time dimension

gdf = ds.vector.to_gdf(reducer=mean)
gdf.head()

	geometry	waterlevel
stations
13670	POINT (12.25342 45.25635)	0.138058
2798	POINT (12.22412 45.27100)	0.140681
2799	POINT (12.29736 45.22705)	0.133489
13775	POINT (12.75879 45.24902)	0.121096
13723	POINT (12.50244 45.25635)	0.126971

# Write the Dataset to an ogr compliant netcdf
from os.path import join

ds.vector.to_netcdf(join("tmpdir", "ds_ogr.nc"), ogr_compliant=True)

# It is indeed ogr compliant as ogrinfo.exe is able to read it
!ogrinfo tmpdir/ds_ogr.nc

Warning 1: The dataset has several variables that could be identified as vector fields, but not all share the same primary dimension. Consequently they will be ignored.
Warning 1: The dataset has several variables that could be identified as vector fields, but not all share the same primary dimension. Consequently they will be ignored.
INFO: Open of `tmpdir/ds_ogr.nc'
      using driver `netCDF' successful.
Metadata:
  NC_GLOBAL#Conventions=CF-1.6
  NC_GLOBAL#coordinates=ogc_wkt spatial_ref
  NC_GLOBAL#GDAL=GDAL 3.8.4
  NC_GLOBAL#ogr_geometry_field=ogc_wkt
  NC_GLOBAL#ogr_layer_type=Point