Projection and placement of time series

[Lucas]

I am new to working with .nc files. I want to calculate the annual electricity production (AEP) for offshore wind turbines around Ireland based on the spatially resolved wind speed time series data from the New European Wind Atlas. However, I encounter the problem that my resulting AEP is not located correctly after I imported it into QGIS for post-processing.

Here are my working steps without aiming for a correct projection and placement :
1. I download the wind data for each month in 2018 (here displayed the URLs for January and December 2018):
https://wps.neweuropeanwindatlas.eu/api/mesoscale-ts/v1/get-data-bbox?southBoundLatitude=50.5&northBoundLatitude=56.1&westBoundLongitude=-11.6&eastBoundLongitude=-5.2&height=150&variable=WS&dt_start=2018-01-01T00:00:00&dt_stop=2018-01-31T23:59:59
...
https://wps.neweuropeanwindatlas.eu/api/mesoscale-ts/v1/get-data-bbox?southBoundLatitude=50.5&northBoundLatitude=56.1&westBoundLongitude=-11.6&eastBoundLongitude=-5.2&height=150&variable=WS&dt_start=2018-12-01T00:00:00&dt_stop=2018-12-31T23:59:59

2. I utilize a Python script with xarray to calculate the AEP with the power curve of a reference turbine:

# Imports
import xarray as xr
import numpy as np
import windkit as wk
import rasterio
import pandas as pd
from rasterio.transform import from_origin

# Constants
rho_air = 1.215  # kg/m³
D_r = 240  # Rotor diameter in meters
P_r = 15e6  # Rated power in Watts
v_cut_in = 3.0  # m/s
v_cut_out = 25.0  # m/s
v_rated = 10.59  # m/s
A = np.pi * (D_r / 2) ** 2  # m²

# Input
# Loop through each month's wind speed data files
for month in range(1, 13):
    file_path = f"R:\...\Wind speed time series\\mesoscale-ts-2018-{month:02d}.nc"
    
    print(f"Loading wind speed dataset for month {month}...")
    data = xr.open_dataset(file_path)  # Einlesen der Daten als xarray.Dataset
    print(f"Dataset for month {month} loaded successfully.")

    # Resample the data to hourly mean
    data = data.resample(time='h').mean()

    # Load power curve
    power_curve_path = r"R:\...\Reference_turbine_2032.csv"
    print("Loading power curve...")
    power_curve = pd.read_csv(power_curve_path)
    print("Power curve loaded successfully.")

    # Interpolate C_p for the wind speeds
    wind_speeds = power_curve['Wind [m/s]'].values
    C_p_values = power_curve['Power Coefficient C_p [-]'].values

    # Interpolate C_p for each time step and location
    C_p_interpolated = np.interp(data['WS'].values.flatten(), wind_speeds, C_p_values, left=0, right=0)
    C_p_interpolated = xr.DataArray(C_p_interpolated.reshape(data['WS'].shape), dims=data['WS'].dims, coords=data['WS'].coords)

    # Calculate the power output based on wind speed conditions
    data_cond_1 = xr.where((data['WS'] >= v_cut_in) & (data['WS'] < v_rated), 0.5 * rho_air * A * C_p_interpolated * data['WS']**3, 0)
    data_cond_2 = xr.where((data['WS'] >= v_rated) & (data['WS'] < v_cut_out), P_r, 0)

    # Total monthly energy production (MEP) for the month
    MEP_month = data_cond_1 + data_cond_2

    # Save the monthly MEP as a separate variable
    globals()[f'MEP_{month:02d}'] = MEP_month

    print(f"Monthly electricity production for month {month} calculated.")

# Sum the MEPs over time to get total monthly production
MEP_01_total = MEP_01.sum(dim='time')
MEP_02_total = MEP_02.sum(dim='time')
MEP_03_total = MEP_03.sum(dim='time')
MEP_04_total = MEP_04.sum(dim='time')
MEP_05_total = MEP_05.sum(dim='time')
MEP_06_total = MEP_06.sum(dim='time')
MEP_07_total = MEP_07.sum(dim='time')
MEP_08_total = MEP_08.sum(dim='time')
MEP_09_total = MEP_09.sum(dim='time')
MEP_10_total = MEP_10.sum(dim='time')
MEP_11_total = MEP_11.sum(dim='time')
MEP_12_total = MEP_12.sum(dim='time')
# Calculate AEP without losses
AEP_total = (MEP_01_total + MEP_02_total + MEP_03_total + MEP_04_total + MEP_05_total + MEP_06_total + MEP_07_total + MEP_08_total + MEP_09_total + MEP_10_total + MEP_11_total + MEP_12_total) / (10**9)

3. I save the AEP as a geotiff:

# Define output file path
output_file_path = r"R:\...\AEP_total.tif"

# Extract coordinates and pixel size
lon = data['XLON'].values  # Assuming XLON is the longitude coordinate
lat = data['XLAT'].values  # Assuming XLAT is the latitude coordinate
west = lon.min()  # Minimum longitude
east = lon.max()  # Maximum longitude
north = lat.max()  # Maximum latitude
south = lat.min()  # Minimum latitude

# Calculate pixel size based on the shape of the data
pixel_size_x = (east - west) / data.sizes['west_east']
pixel_size_y = (north - south) / data.sizes['south_north']

# Get the transform for the output raster
transform = from_origin(west, north, pixel_size_x, pixel_size_y)

# Squeeze AEP_total to remove extra dimensions
AEP_total_squeezed = AEP_total.squeeze()

# Save the AEP_total as a GeoTIFF
with rasterio.open(
    output_file_path,
    'w',
    driver='GTiff',
    height=AEP_total_squeezed.shape[0],
    width=AEP_total_squeezed.shape[1],
    count=1,
    dtype=AEP_total_squeezed.dtype,
    crs='EPSG:4326',  # Set the CRS based on your requirements
    transform=transform
) as dst:
    dst.write(AEP_total_squeezed.values, 1)

print(f"AEP_total saved to {output_file_path}.")

4. Import the .tif file into QGIS.

Overall i tried several ways to get a correct projection, placement and size of my output data. In the end I need the AEP within EPSG:2157 for post-processing. For example I tired utilizing the windkit package, as well as warping the AEP file with gdal.
However, I was not able to come up with a solution. The output was either located somewhere else in QGIS or was rotated and bigger.

Where in my working process is the correct place to interfere and ensure that the data is processed in the correct way.
Thanks a lot and best regards
Lucas

Edited January 28 by Lucas

[Rogier]

The coordinates from the original WRF projection are already stored in the west_east and south_north coordinates and the spatial reference is already stored in the crs variable. So to make a tif you can do:

import xarray as xr

ds = xr.open_dataset("mesoscale-ts.nc")

ds_proj_mean = ds.WS.mean("time")

ds_proj_mean.rename({"crs":"spatial_ref"}).rio.to_raster("ws_150.tif")

That opens it in the right place for me:

or to get it in 2157:

ds_proj_mean = ds.WS.mean("time")

ds_proj_mean = ds_proj_mean.rename({"crs":"spatial_ref"}).drop_vars(["XLAT","XLON"]).rio.reproject("EPSG:2157")

ds_proj_mean.rio.to_raster("ws_150_2157.tif")

Edited February 3 by Rogier

[Lucas]

Thanks Rogier, that resolves the issue.