Geometries

See also

1. The CF conventions on Geometries

2. cf_xarray.shapely_to_cf()

3. cf_xarray.cf_to_shapely()

cf_xarray can convert between vector geometries represented as shapely objects and CF-compliant array representations of those geometries.

Let’s start by creating an xarray object containing some shapely geometries. This example uses a xr.DataArray but these functions also work with a xr.Dataset where one of the data variables contains an array of shapes.

import cf_xarray as cfxr
import xarray as xr

from shapely.geometry import MultiPoint, Point

da = xr.DataArray(
    [
        MultiPoint([(1.0, 2.0), (2.0, 3.0)]),
        Point(3.0, 4.0),
        Point(4.0, 5.0),
        Point(3.0, 4.0),
    ],
    dims=("index",),
    name="geometry"
)

Warning

cf_xarray does not support handle multiple types of shapes (Point, Line, Polygon) in one xr.DataArray, but multipart geometries are supported and can be mixed with single-part geometries of the same type.

Now we can take that xr.DataArray containing shapely geometries and convert it to cf:

ds_cf = cfxr.shapely_to_cf(da)
ds_cf

<xarray.Dataset> Size: 184B
Dimensions:             (index: 4, node: 5)
Coordinates:
    x                   (node) float64 40B 1.0 2.0 3.0 4.0 3.0
    y                   (node) float64 40B 2.0 3.0 4.0 5.0 4.0
    crd_x               (index) float64 32B 1.0 3.0 4.0 3.0
    crd_y               (index) float64 32B 2.0 4.0 5.0 4.0
Dimensions without coordinates: index, node
Data variables:
    node_count          (index) int64 32B 2 1 1 1
    geometry_container  float64 8B nan

xarray.Dataset

• Dimensions:
  • index: 4
  • node: 5
• Coordinates: (4)
  • x
    (node)
    float64
    1.0 2.0 3.0 4.0 3.0

    axis :

    X

    array([1., 2., 3., 4., 3.])

  • y
    (node)
    float64
    2.0 3.0 4.0 5.0 4.0

    axis :

    Y

    array([2., 3., 4., 5., 4.])

  • crd_x
    (index)
    float64
    1.0 3.0 4.0 3.0

    nodes :

    x

    array([1., 3., 4., 3.])

  • crd_y
    (index)
    float64
    2.0 4.0 5.0 4.0

    nodes :

    y

    array([2., 4., 5., 4.])

• Data variables: (2)
  • node_count
    (index)
    int64
    2 1 1 1

    array([2, 1, 1, 1])

  • geometry_container
    ()
    float64
    nan

    geometry_type :

    point

    node_count :

    node_count

    node_coordinates :

    x y

    coordinates :

    crd_x crd_y

    array(nan)

• Indexes: (0)
• Attributes: (0)

This function returns a xr.Dataset containing the CF fields needed to reconstruct the geometries. In particular there are:

• 'x', 'y' : the node coordinates

• 'crd_x', 'crd_y' : the feature coordinates (might have different names if grid_mapping is available).

• 'node_count' : The number of nodes per feature. Always present for Lines and Polygons. For Points: only present if there are multipart geometries.

• 'part_node_count' : The number of nodes per individual geometry. Only for Lines with multipart geometries and for Polygons with multipart geometries or holes.

• 'interior_ring' : Integer boolean indicating whether ring is interior or exterior. Only for Polygons with holes.

• 'geometry_container’ : Empty variable with attributes describing the geometry type.

Here are the attributes on geometry_container. This pattern mimics the convention of specifying spatial reference information in the attrs of the empty array spatial_ref.

ds_cf.geometry_container.attrs

{'geometry_type': 'point',
 'node_count': 'node_count',
 'node_coordinates': 'x y',
 'coordinates': 'crd_x crd_y'}

Note

Z axis is not yet supported for any shapes.

This xr.Dataset can be converted back into a xr.DataArray of shapely geometries:

cfxr.cf_to_shapely(ds_cf)

<xarray.DataArray 'geometry' (index: 4)> Size: 32B
array([<MULTIPOINT ((1 2), (2 3))>, <POINT (3 4)>, <POINT (4 5)>,
       <POINT (3 4)>], dtype=object)
Coordinates:
    crd_x    (index) float64 32B 1.0 3.0 4.0 3.0
    crd_y    (index) float64 32B 2.0 4.0 5.0 4.0
Dimensions without coordinates: index

xarray.DataArray

'geometry'

• index: 4

• MULTIPOINT ((1 2), (2 3)) POINT (3 4) POINT (4 5) POINT (3 4)

  array([<MULTIPOINT ((1 2), (2 3))>, <POINT (3 4)>, <POINT (4 5)>,
         <POINT (3 4)>], dtype=object)

• Coordinates: (2)
  • crd_x
    (index)
    float64
    1.0 3.0 4.0 3.0

    nodes :

    x

    array([1., 3., 4., 3.])

  • crd_y
    (index)
    float64
    2.0 4.0 5.0 4.0

    nodes :

    y

    array([2., 4., 5., 4.])

• Indexes: (0)
• Attributes: (0)

This conversion adds coordinates that aren’t in the xr.DataArray that we started with. By default these are called crd_x and crd_y unless grid_mapping is specified.

Gotchas

For MultiPolygons with holes the CF notation is slightly ambiguous on which hole is associated with which polygon. This is problematic because shapely stores holes within the polygon object that they are associated with. cf_xarray assumes that the the shapes are interleaved such that the holes (interior rings) are associated with the exteriors (exterior rings) that immediately precede them.