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 from UTM zone 32N to WGS84 with bilinear interpolation
warped = wk.spatial.warp(raster_ds, to_crs="EPSG:4326", method="bilinear")
print(f"Warped shape: {dict(warped.sizes)}")
print(f"Longitude range: {warped.west_east.min().item():.3f} to {warped.west_east.max().item():.3f}")

Warped shape: {'height': 1, 'south_north': 70, 'west_east': 125}
Longitude range: 9.000 to 9.161

The method is passed on to rasterio.warp.reproject(), so any of its resampling methods ("nearest", "bilinear", "cubic", …) can be used.

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=502_000, miny=6_202_000,
    maxx=508_000, maxy=6_208_000,
    crs="EPSG:32632"
)

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

Original shape: {'height': 1, 'south_north': 101, 'west_east': 101}
Clipped shape: {'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(504_000, 506_000, 21),
    south_north=np.linspace(6_204_000, 6_206_000, 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: {dict(clipped_margin.sizes)}")

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

Masking

windkit.spatial.mask()

Set values outside a geometry to NaN, keeping the original shape (unlike clip, which drops them):

# Mask the raster with the bounding box from above
masked = wk.spatial.mask(raster_ds, bbox)
print(f"Masked shape: {dict(masked.sizes)}")
print(f"Cells kept: {masked.output.notnull().sum().item()} of {masked.output.size}")

Masked shape: {'height': 1, 'south_north': 101, 'west_east': 101}
Cells kept: 3721 of 10201

The mask can also be a geopandas.GeoDataFrame, for example a site boundary polygon. Pass invert=True to null the values inside the mask instead.

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, natural, rbf
)

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"
)

Radial basis function interpolation

The linear and cubic methods triangulate the source points and return NaN outside their convex hull, while nearest takes the value of the closest source point everywhere. method="rbf" instead fits a radial basis function through the source points with scipy.interpolate.RBFInterpolator, which differs in three ways:

• It extrapolates smoothly outside the convex hull instead of returning NaN.

• It works with very few source points. Two points can be interpolated with kernel="linear"; the default thin_plate_spline kernel needs at least three points that are not collinear, because its polynomial term is of degree 1. Collinear sources, such as masts along a transect, need kernel="linear" or degree=0.

• It interpolates horizontally only. Radial distances are measured in the coordinate units, so pass exclude_dims=["height"] and use a projected CRS. A geographic CRS raises a warning because degrees of longitude and latitude are not the same length.

height can only be excluded when it is a dimension of the source, which is the case for stacked_point and cuboid structures. In point structure every (location, height) pair is its own point, so a source with several heights per location has duplicate horizontal coordinates and is rejected. Convert such a source with windkit.spatial.to_stacked_point() first.

Because the fit is global, a NaN value at one source point makes the whole slice it belongs to (for example one sector) NaN at every target. Drop or fill missing values before interpolating.

Extra keyword arguments are passed on to RBFInterpolator: kernel, smoothing, epsilon, degree and neighbors.

# Interpolate between a few met masts with a linear kernel
interpolated = wk.spatial.interp_unstructured_like(
    wk.spatial.to_stacked_point(mast_data),
    target_points,
    exclude_dims=["height"],
    method="rbf",
    kernel="linear",
)

# Smooth a large point cloud with a local fit
interpolated = wk.spatial.interp_unstructured_like(
    point_cloud,
    target_grid,
    exclude_dims=["height"],
    method="rbf",
    smoothing=1.0,
    neighbors=50,
)

Note

The default global fit solves a dense linear system whose cost grows with the cube of the number of source points, and becomes slow beyond a few thousand points. Pass neighbors to fit locally in that case.

Nearest Points

windkit.spatial.nearest_points()

Select the nearest source cell for every target location without interpolating source values:

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",
)
source["wind_speed"] = ("point", [7.0, 8.0, 9.0])
target = wk.spatial.create_point(
    west_east=[130.0, 275.0],
    south_north=[130.0, 275.0],
    height=[80.0, 80.0],
    crs="EPSG:32632",
)

selected = wk.spatial.nearest_points(
    source, target, return_distance=True
)
print(selected[["wind_speed", "distance"]])

<xarray.Dataset> Size: 81B
Dimensions:      (point: 2)
Coordinates:
    height       (point) float64 16B 80.0 80.0
    south_north  (point) float64 16B 100.0 300.0
    west_east    (point) float64 16B 100.0 300.0
    crs          int8 1B 0
Dimensions without coordinates: point
Data variables:
    wind_speed   (point) float64 16B 7.0 9.0
    distance     (point) float64 16B 42.43 35.36
Attributes:
    Conventions:  CF-1.8
    history:      2026-09-29T09:46:04+00:00:\twindkit==2.2.0\t"point")\n2026-...

The result contains exact coordinates and values selected from source, using the spatial structure of target. If the exact source coordinates cannot form that structure, nearest_points raises ValueError rather than changing the return type. Convert the target explicitly with wk.spatial.to_point(target) to request lossless point output. Lookup dimensions can be selected with dims. Omitting a source spatial coordinate selects its first representative and emits a warning.

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-09-29T09:46:04+00:00:\twindkit==2.2.0\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-09-29T09:46:04+00:00: windkit==2.2.0 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)