Operations

WindKit provides functions for common geospatial operations on xarray datasets.

Reprojection

Transform data between coordinate reference systems.

windkit.spatial.reproject()

Reproject point-based datasets to a new CRS:

# Reproject point data from WGS84 to UTM zone 32N
point_utm = wk.spatial.reproject(point_ds, to_crs="EPSG:32632")
print(f"Original CRS: {wk.spatial.get_crs(point_ds)}")
print(f"New CRS: {wk.spatial.get_crs(point_utm)}")
print(f"Original coords: west_east={point_ds.west_east.values}")
print(f"Reprojected coords: west_east={point_utm.west_east.values}")

Original CRS: 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)"],MEMBER["World Geodetic System 1984 (G2296)"],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]]
New CRS: PROJCRS["WGS 84 / UTM zone 32N",BASEGEOGCRS["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)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 32N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State."],BBOX[0,6,84,12]],ID["EPSG",32632]]
Original coords: west_east=[10. 11.]
Reprojected coords: west_east=[562366.63018029 621485.72767786]

windkit.spatial.warp()

Reproject and resample raster data to a new CRS:

# Warp raster to new CRS with bilinear interpolation
# warped = wk.spatial.warp(raster_ds, to_crs="EPSG:4326", method="bilinear")
# Available methods: nearest, bilinear, cubic, lanczos
print("Raster warping: wk.spatial.warp(raster, to_crs='EPSG:4326', method='bilinear')")

Raster warping: wk.spatial.warp(raster, to_crs='EPSG:4326', method='bilinear')

Clipping

Extract data within a spatial extent.

windkit.spatial.clip()

Clip data to a bounding box:

# Create a bounding box
bbox = wk.spatial.BBox.from_cornerpts(
    minx=2000, miny=2000,
    maxx=8000, maxy=8000,
    crs="EPSG:32632"
)

# Clip raster to bounding box
clipped = wk.spatial.clip(raster_ds, bbox)
print(f"Original shape: {raster_ds.dims}")
print(f"Clipped shape: {clipped.dims}")

Original shape: FrozenMappingWarningOnValuesAccess({'height': 1, 'south_north': 101, 'west_east': 101})
Clipped shape: FrozenMappingWarningOnValuesAccess({'height': 1, 'south_north': 61, 'west_east': 61})

windkit.spatial.clip_with_margin()

Clip with a margin around another dataset (margin is a multiple of grid spacing):

# Create a smaller target dataset
target_raster = wk.spatial.create_raster(
    west_east=np.linspace(4000, 6000, 21),
    south_north=np.linspace(4000, 6000, 21),
    crs="EPSG:32632"
)

# Clip with margin (default 5x grid spacing around target)
clipped_margin = wk.spatial.clip_with_margin(
    raster_ds,
    target_raster,
    margin_dx_factor=3.0
)
print(f"Clipped with margin shape: {clipped_margin.dims}")

Clipped with margin shape: FrozenMappingWarningOnValuesAccess({'height': 1, 'south_north': 27, 'west_east': 27})

Masking

windkit.spatial.mask()

Apply a spatial mask to data:

# Masking example - mask data outside a geometry
# masked = wk.spatial.mask(dataset, geometry)
print("Spatial masking: wk.spatial.mask(dataset, geometry)")

Spatial masking: wk.spatial.mask(dataset, geometry)

Interpolation

Interpolate data to new spatial locations. These functions require datasets with actual data variables (not just coordinates).

windkit.spatial.interp_structured_like()

Interpolate structured (raster/cuboid) data to match another dataset’s grid:

# Interpolate source data to match target grid resolution
interpolated = wk.spatial.interp_structured_like(
    source_data,
    target_grid,
    method="linear"  # Options: nearest, linear
)

windkit.spatial.interp_unstructured()

Interpolate structured data to unstructured (point) locations:

import xarray as xr

# Create target coordinates as DataArrays (must share a common dimension)
we = xr.DataArray([2500.0, 5000.0, 7500.0], dims="point")
sn = xr.DataArray([2500.0, 5000.0, 7500.0], dims="point")

# Interpolate to specific point locations
interpolated = wk.spatial.interp_unstructured(
    source_data,
    west_east=we,
    south_north=sn,
    method="linear"  # Options: nearest, linear, cubic
)

windkit.spatial.interp_unstructured_like()

Interpolate to match unstructured point locations from another dataset:

# Interpolate source data to match target point locations
interpolated = wk.spatial.interp_unstructured_like(
    source_data,
    target_points,
    method="linear"
)

Nearest Points

windkit.spatial.nearest_points()

Find nearest points between two datasets:

# Create two point datasets
source = wk.spatial.create_point(
    west_east=[100.0, 200.0, 300.0],
    south_north=[100.0, 200.0, 300.0],
    height=[80.0, 80.0, 80.0],
    crs="EPSG:32632"
)
target = wk.spatial.create_point(
    west_east=[150.0, 250.0],
    south_north=[150.0, 250.0],
    height=[80.0, 80.0],
    crs="EPSG:32632"
)

# Find nearest points in target for each point in source
indices = wk.spatial.nearest_points(source, target)
print(f"Nearest target indices for each source point: {indices}")

Nearest target indices for each source point: <xarray.Dataset> Size: 65B
Dimensions:      (point: 2)
Coordinates:
    height       (point) float64 16B 80.0 80.0
    south_north  (point) float64 16B 100.0 200.0
    west_east    (point) float64 16B 100.0 200.0
    crs          int8 1B 0
Dimensions without coordinates: point
Data variables:
    output       (point) float64 16B 0.0 0.0
Attributes:
    Conventions:  CF-1.8
    history:      2026-01-29T09:29:59+00:00:\twindkit==2.0.1\t"point")\n2026-...

Spatial Comparison

windkit.spatial.are_spatially_equal()

Check if two datasets have identical spatial coordinates:

ds1 = wk.spatial.create_point(
    west_east=[10.0], south_north=[56.0], height=[80.0], crs="EPSG:4326"
)
ds2 = ds1.copy()
ds3 = wk.spatial.create_point(
    west_east=[11.0], south_north=[57.0], height=[80.0], crs="EPSG:4326"
)

print(f"ds1 == ds2: {wk.spatial.are_spatially_equal(ds1, ds2)}")
print(f"ds1 == ds3: {wk.spatial.are_spatially_equal(ds1, ds3)}")

ds1 == ds2: True
ds1 == ds3: False

windkit.spatial.equal_spatial_shape()

Check if two datasets have the same spatial shape:

raster1 = wk.spatial.create_raster(
    west_east=np.linspace(0, 100, 11),
    south_north=np.linspace(0, 100, 11),
    crs="EPSG:32632"
)
raster2 = wk.spatial.create_raster(
    west_east=np.linspace(0, 200, 11),  # Different extent, same shape
    south_north=np.linspace(0, 200, 11),
    crs="EPSG:32632"
)

same_shape = wk.spatial.equal_spatial_shape(raster1, raster2)
print(f"Same spatial shape: {same_shape}")

Same spatial shape: True

windkit.spatial.covers()

Check if one dataset spatially covers another:

large = wk.spatial.create_raster(
    west_east=np.linspace(0, 1000, 101),
    south_north=np.linspace(0, 1000, 101),
    crs="EPSG:32632"
)
small = wk.spatial.create_raster(
    west_east=np.linspace(200, 800, 61),
    south_north=np.linspace(200, 800, 61),
    crs="EPSG:32632"
)

print(f"Large covers small: {wk.spatial.covers(large, small)}")
print(f"Small covers large: {wk.spatial.covers(small, large)}")

Large covers small: <xarray.DataArray ()> Size: 1B
array(True)
Coordinates:
    crs      int8 1B 0
Attributes:
    URI:            
    units:          m
    axis:           Z
    description:    It is the height above ground either related to the measu...
    long_name:      Height above ground
    standard_name:  height
Small covers large: <xarray.DataArray ()> Size: 1B
array(False)
Coordinates:
    crs      int8 1B 0
Attributes:
    URI:            
    units:          m
    axis:           Z
    description:    It is the height above ground either related to the measu...
    long_name:      Height above ground
    standard_name:  height

Counting Points

windkit.spatial.count_spatial_points()

Count the number of spatial points in a dataset:

n_points = wk.spatial.count_spatial_points(point_ds)
print(f"Number of spatial points: {n_points}")

Number of spatial points: 2

GeoDataFrame Conversion

Convert between xarray datasets and geopandas GeoDataFrames.

windkit.spatial.gdf_to_ds()

Convert GeoDataFrame to xarray Dataset:

import geopandas as gpd
from shapely.geometry import Point

# Create a GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"value": [1, 2, 3]},
    geometry=[Point(10, 56), Point(11, 57), Point(12, 58)],
    crs="EPSG:4326"
)

# Convert to xarray Dataset
ds_from_gdf = wk.spatial.gdf_to_ds(gdf, height=80.0)
ds_from_gdf

<xarray.Dataset> Size: 129B
Dimensions:      (height: 1, south_north: 3, west_east: 3)
Coordinates:
  * height       (height) float64 8B 80.0
  * south_north  (south_north) float64 24B 56.0 57.0 58.0
  * west_east    (west_east) float64 24B 10.0 11.0 12.0
    crs          int8 1B 0
Data variables:
    output       (height, south_north, west_east) float64 72B 0.0 0.0 ... 0.0
Attributes:
    Conventions:  CF-1.8
    history:      2026-01-29T09:29:59+00:00:\twindkit==2.0.1\tcreate_dataset(...

xarray.Dataset

• Dimensions:
  • height: 1
  • south_north: 3
  • west_east: 3
• Coordinates: (4)
  • height
    (height)
    float64
    80.0

    URI :

    axis :

    Z

    description :

    It is the height above ground either related to the measurements or to the hub of a turbine.

    long_name :

    Height above ground

    standard_name :

    height

    units :

    m

    array([80.])

  • south_north
    (south_north)
    float64
    56.0 57.0 58.0

    URI :

    axis :

    Y

    description :

    Positive northward direction in geographic coordinates (degrees).

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([56., 57., 58.])

  • west_east
    (west_east)
    float64
    10.0 11.0 12.0

    URI :

    axis :

    X

    description :

    Positive eastward direction in geographic coordinates (degrees).

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([10., 11., 12.])

  • crs
    ()
    int8
    0

    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)"],MEMBER["World Geodetic System 1984 (G2296)"],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

    long_name :

    CRS definition

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

    array(0, dtype=int8)

• Data variables: (1)
  • output
    (height, south_north, west_east)
    float64
    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

    array([[[0., 0., 0.],
            [0., 0., 0.],
            [0., 0., 0.]]])

• Attributes: (2)

  Conventions :

  CF-1.8

  history :

  2026-01-29T09:29:59+00:00: windkit==2.0.1 create_dataset(x, y, z, struct=struct, crs=gdf.crs)

windkit.spatial.ds_to_gdf()

Convert xarray Dataset to GeoDataFrame:

# Convert back to GeoDataFrame
gdf_from_ds = wk.spatial.ds_to_gdf(point_ds)
gdf_from_ds

	geometry
0	POINT (10 56)
1	POINT (11 57)