Add minimal example

example/README.md

@@ -0,0 +1,29 @@
+# Minimal example for computing capacity factors using atlite
+
+This is a minimal example that produces capacity factors for pv and wind
+for the Netherlands. The input data has been cropped and downsampled to keep the files small.
+
+To run the example, first create the conda environment
+
+    conda env create -f environment.yaml
+
+To download the ERA5 cutout, you need to provide your credentials in a file
+called .cdsapirc. Given that, you can run
+
+    python example/download_cutout_era5.py
+
+, which will save the cutout in example/results/cutout_era5.nc. Next, you can run
+
+    prepare_cf_from_raster_layout.py
+
+and will find two files, example/results/cf_pv.nc and example/results/cf_wind.nc, containing the capacity factors.
+
+## Use cases
+
+Provide a raster layout, aggregate results to polygons. See `prepare_cf_from_raster_layout.py`
+
+Provide a point layout, aggregate results to polygons. See `prepare_cf_from_point_layout.py`
+
+Provide no layout, get capacity factors as raster at the resolution of the cutout. Same as one of the above, but do not pass `shapes`, `layout`, `matrix`, `per_unit`. Instead, pass `capacity_factor_timeseries=True`. See `prepare_cf_from_raster_no_layout_to_raster.py`.
+
+Provide a raster or point layout, do not aggregate, get capacity factors as raster at cutout resolution. See `prepare_cf_from_raster_layout_to_raster.py`.

example/data/layout_points_NLD.csv

@@ -0,0 +1,3 @@
+id,techs,lat,lon,capacity
+0,wind,51.956667748480086,4.320588199650519,2
+1,wind,51.98127968603208,4.405069363584079,4

example/download_cutout_era5.py

@@ -0,0 +1,20 @@
+import atlite
+import logging
+from pathlib import Path
+
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    here = Path(__file__).parent
+    path_cutout = here / "results" / "cutout_era5.nc"
+
+    features = None
+    cutout = atlite.Cutout(
+        path_cutout, 
+        module="era5",
+        x=slice(3, 7.5),
+        y=slice(50.5, 54),
+        time=slice("2013-01-01", "2013-01-02"),
+    )
+    cutout.prepare(features=features)

example/environment.yaml

@@ -0,0 +1,10 @@
+name: module_pv_wind_example
+channels:
+    - conda-forge
+    - bioconda    
+dependencies:
+    - rioxarray
+    - matplotlib
+    - geopandas
+    - atlite
+    - cdsapi

example/plot_cf.py

@@ -0,0 +1,48 @@
+#%% 
+import xarray as xr
+import pandas as pd
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+# load data
+here = Path(__file__).parent
+
+cf_wind_layout_point = xr.open_dataset(here / "results" / "cf_wind_layout_point.nc")
+
+cf_wind_layout_raster = xr.open_dataset(here / "results" / "cf_wind_layout_raster.nc")
+
+cf_wind_layout_raster_to_raster = xr.open_dataset(here / "results" / "cf_wind_layout_raster_to_raster.nc")
+cf_wind_layout_raster_to_raster
+
+#%%
+df_cf = cf_wind_layout_point.to_dataframe()
+
+idx=pd.IndexSlice
+
+fig,ax=plt.subplots()
+for id in df_cf.index.get_level_values("id").unique():
+    timeseries = df_cf.loc[idx[:, id], "__xarray_dataarray_variable__"]
+    timeseries.plot(ax=ax, label=id)
+
+plt.legend(
+    loc="upper right",
+    bbox_to_anchor=(1.2, 1),
+    bbox_transform=ax.transAxes,
+)
+
+#%%
+df_cf = cf_wind_layout_raster.to_dataframe()
+
+fig,ax=plt.subplots()
+for id in df_cf.index.get_level_values("id").unique():
+    timeseries = df_cf.loc[idx[:, id], "__xarray_dataarray_variable__"]
+    timeseries.plot(ax=ax, label=id)
+
+plt.legend(
+    loc="upper right",
+    bbox_to_anchor=(1.2, 1),
+    bbox_transform=ax.transAxes,
+)
+
+#%%
+cf_wind_layout_raster_to_raster

example/prepare_cf_from_point_layout.py

@@ -0,0 +1,48 @@
+import atlite
+import rioxarray as rxr
+import rasterio as rio
+import pandas as pd
+import geopandas as gpd
+from pathlib import Path
+
+
+def read_yaml(filepath):
+    with open(filepath, "r") as file:
+        return yaml.safe_load(file)
+
+
+if __name__ == "__main__":
+    here = Path(__file__).parent
+    path_layout = here / "data" / "layout_points_NLD.csv"
+    path_spatial_units = here / "data" / "regional_NLD.geojson"
+    path_cutout = here / "results" / "cutout_era5.nc"
+    path_cf_pv = here / "results" / "cf_pv_layout_point.nc"
+    path_cf_wind = here / "results" / "cf_wind_layout_point.nc"
+
+    # load data
+    cutout = atlite.Cutout(path_cutout)
+    layout_point = pd.read_csv(path_layout)
+    spatial_units = gpd.read_file(path_spatial_units).set_index("id")
+
+
+    # prepare layout from list of points
+    layout_point = layout_point.rename(columns={"lon": "x", "lat": "y"})
+    layout = cutout.layout_from_capacity_list(layout_point, col="capacity")
+
+    # compute capacity factors
+    capacityfactors_pv = cutout.pv(
+        panel="CSi", 
+        orientation={"slope": 30.0, "azimuth": 180.0}, 
+        layout=layout,
+        shapes=spatial_units,
+        per_unit=True,
+    )
+    capacityfactors_pv.to_netcdf(path_cf_pv)
+
+    capacityfactors_wind = cutout.wind(
+        turbine="Vestas_V90_3MW", 
+        layout=layout,
+        shapes=spatial_units,
+        per_unit=True,
+    )
+    capacityfactors_wind.to_netcdf(path_cf_wind)

example/prepare_cf_from_raster_layout.py

@@ -0,0 +1,51 @@
+import atlite
+import rioxarray as rxr
+import rasterio as rio
+import geopandas as gpd
+from pathlib import Path
+
+
+def read_yaml(filepath):
+    with open(filepath, "r") as file:
+        return yaml.safe_load(file)
+
+
+if __name__ == "__main__":
+    here = Path(__file__).parent
+    path_availability = here / "data" / "availability_NLD.tif"
+    path_spatial_units = here / "data" / "regional_NLD.geojson"
+    path_cutout = here / "results" / "cutout_era5.nc"
+    path_cf_pv = here / "results" / "cf_pv_layout_raster.nc"
+    path_cf_wind = here / "results" / "cf_wind_layout_raster.nc"
+
+    # load data
+    cutout = atlite.Cutout(path_cutout)
+    availability = rxr.open_rasterio(path_availability)
+    spatial_units = gpd.read_file(path_spatial_units).set_index("id")
+
+    # resample availability to the resolution of the cutout
+    match = cutout.uniform_layout().rio.write_crs(cutout.crs).rio.write_transform(cutout.transform)
+    layout = availability.squeeze(drop=True)
+    layout = layout.rio.reproject_match(
+        match,
+        resampling=rio.enums.Resampling.sum,
+        nodata=0
+    )
+
+    # compute capacity factors
+    capacityfactors_pv = cutout.pv(
+        panel="CSi", 
+        orientation={"slope": 30.0, "azimuth": 180.0}, 
+        layout=layout,
+        shapes=spatial_units,
+        per_unit=True,
+    )
+    capacityfactors_pv.to_netcdf(path_cf_pv)
+
+    capacityfactors_wind = cutout.wind(
+        turbine="Vestas_V90_3MW", 
+        layout=layout,
+        shapes=spatial_units,
+        per_unit=True,
+    )
+    capacityfactors_wind.to_netcdf(path_cf_wind)

example/prepare_cf_from_raster_layout_to_raster.py

@@ -0,0 +1,55 @@
+import atlite
+import rioxarray as rxr
+import rasterio as rio
+import geopandas as gpd
+from pathlib import Path
+
+
+def read_yaml(filepath):
+    with open(filepath, "r") as file:
+        return yaml.safe_load(file)
+
+
+if __name__ == "__main__":
+    here = Path(__file__).parent
+    path_availability = here / "data" / "availability_NLD.tif"
+    path_spatial_units = here / "data" / "regional_NLD.geojson"
+    path_cutout = here / "results" / "cutout_era5.nc"
+    path_cf_pv = here / "results" / "cf_pv_layout_raster_to_raster.nc"
+    path_cf_wind = here / "results" / "cf_wind_layout_raster_to_raster.nc"
+
+    # load data
+    cutout = atlite.Cutout(path_cutout)
+    availability = rxr.open_rasterio(path_availability)
+    spatial_units = gpd.read_file(path_spatial_units).set_index("id")
+
+    # resample availability to the resolution of the cutout
+    match = cutout.uniform_layout().rio.write_crs(cutout.crs).rio.write_transform(cutout.transform)
+    layout = availability.squeeze(drop=True)
+    layout = layout.rio.reproject_match(
+        match,
+        resampling=rio.enums.Resampling.sum,
+        nodata=0
+    )
+
+    # experimental/hacky: create a matrix with an extra dimension with a coordinate 
+    # for each pixel to "aggregate" to each pixel, which means not aggregating.
+    matrix = layout.stack(spatial=["y", "x"])
+    matrix = matrix.expand_dims(dim={"name":matrix.spatial.values}, axis=0)
+    matrix = matrix.data
+
+    # compute capacity factors
+    capacityfactors_pv = cutout.pv(
+        panel="CSi", 
+        orientation={"slope": 30.0, "azimuth": 180.0}, 
+        matrix=matrix,
+        per_unit=True,
+    )
+    capacityfactors_pv.to_netcdf(path_cf_pv)
+
+    capacityfactors_wind = cutout.wind(
+        turbine="Vestas_V90_3MW", 
+        matrix=matrix,
+        per_unit=True,
+    )
+    capacityfactors_wind.to_netcdf(path_cf_wind)

example/prepare_cf_from_raster_no_layout_to_raster.py

@@ -0,0 +1,53 @@
+import atlite
+import rioxarray as rxr
+import rasterio as rio
+import geopandas as gpd
+from pathlib import Path
+
+
+def read_yaml(filepath):
+    with open(filepath, "r") as file:
+        return yaml.safe_load(file)
+
+
+if __name__ == "__main__":
+    here = Path(__file__).parent
+    path_availability = here / "data" / "availability_NLD.tif"
+    path_spatial_units = here / "data" / "regional_NLD.geojson"
+    path_cutout = here / "data" / "era5.nc"
+    path_cf_pv = here / "results" / "cf_pv_no_layout_raster_to_raster.nc"
+    path_cf_wind = here / "results" / "cf_wind_no_layout_raster_to_raster.nc"
+
+    # load data
+    cutout = atlite.Cutout(path_cutout)
+    availability = rxr.open_rasterio(path_availability)
+    spatial_units = gpd.read_file(path_spatial_units).set_index("id")
+
+    # resample availability to the resolution of the cutout
+    match = cutout.uniform_layout().rio.write_crs(cutout.crs).rio.write_transform(cutout.transform)
+    layout = availability.squeeze(drop=True)
+    layout = layout.rio.reproject_match(
+        match,
+        resampling=rio.enums.Resampling.sum,
+        nodata=0
+    )
+
+    # experimental/hacky: create a matrix with an extra dimension with a coordinate 
+    # for each pixel to "aggregate" to each pixel, which means not aggregating.
+    matrix = layout.stack(spatial=["y", "x"])
+    matrix = matrix.expand_dims(dim={"name":matrix.spatial.values}, axis=0)
+    matrix = matrix.data
+
+    # compute capacity factors
+    capacityfactors_pv = cutout.pv(
+        panel="CSi", 
+        orientation={"slope": 30.0, "azimuth": 180.0}, 
+        capacity_factor_timeseries=True
+    )
+    capacityfactors_pv.to_netcdf(path_cf_pv)
+
+    capacityfactors_wind = cutout.wind(
+        turbine="Vestas_V90_3MW", 
+        capacity_factor_timeseries=True
+    )
+    capacityfactors_wind.to_netcdf(path_cf_wind)