Handle Rotation Inside Aeolis to Ensure Correct Georeferencing

[frederikvand]

Description

Currently, geospatial data used as input for aeolis-python sometimes has missing data due to a non-rectangular spatial extent during land survey sampling. To address this, Aeolis users often preprocess arrays by rotating and resampling raster data so they align with the coastline to ensure that there is no missing data. However, this potential rotation is currently not stored in the output netCDF file, and arrays are not rotated back, which leads to faulty georeferencing when the netCDF output is interpreted or linked to other models.

Solution Overview

To handle the rotation of input rasters within aeolis-python, I propose the following enhancements:

1. Add a configuration file parameter for specifying the rotation angle of the input raster. This will allow users to provide the rotation information explicitly. The rotation parameter should be added to the aeolis/constants.py file.

2. Integrate functions to handle the rotation of the input raster based on the specified rotation angle. These functions should be added to the aeolis/utils.py file. This will ensure that the model uses the rotated raster internally while preserving the original orientation information.

3. Rotate the output arrays back to the original orientation before writing them to the netCDF file. This will ensure proper georeferencing and compatibility with other models and geospatial packages.

4. Store the original geotiff metadata in the output netCDF file. This will maintain the traceability of the input data and enable proper interpretation of the model results. This enhancement is related to the issue "Integrate Reading Geotiffs for Improved Remote Sensing Data Integration", where I discuss improving the integration of geotiff data.

It is important to note that I do not propose geospatial rotation functions here provided by libraries like rasterio or GDAL. "Rasterio transform" rotates the geospatial metadata but does not change the underlying array data. In the case of aeolis, it is required for the input array itself to be rotated, and the missing data needs to be removed to ensure proper functioning of the model. Therefore, I propose to implement custom rotation functions that operate directly on the array data.

Required Functions

The following functions will be needed to implement the proposed solution and should be added to the aeolis/utils.py file:

import numpy as np
from scipy.ndimage import rotate

def rotate_array(array, angle):
    """
    Rotates the given array by the specified angle.

    Parameters:
    - array (numpy.ndarray): The input array to be rotated.
    - angle (float): The rotation angle in degrees.

    Returns:
    - rotated_array (numpy.ndarray): The rotated array.
    """
    rotated_array = rotate(array, angle, reshape=False, order=1)
    return rotated_array

def remove_missing_data(array):
    """
    Removes rows containing missing data (NaN values) from the array and stores the indices of the removed rows.

    Parameters:
    - array (numpy.ndarray): The input array.

    Returns:
    - cropped_array (numpy.ndarray): The array with missing data removed.
    - removed_indices (numpy.ndarray): The indices of the removed rows.
    """
    mask = np.isnan(array)
    valid_rows = np.all(~mask, axis=1)
    removed_indices = np.where(~valid_rows)[0]
    cropped_array = array[valid_rows]
    return cropped_array, removed_indices

def add_missing_rows(array, removed_indices, original_shape):
    """
    Adds back the removed rows as NaN values to match the original shape of the array.

    Parameters:
    - array (numpy.ndarray): The input array.
    - removed_indices (numpy.ndarray): The indices of the removed rows.
    - original_shape (tuple): The original shape of the array before removing missing data.

    Returns:
    - padded_array (numpy.ndarray): The array with missing rows added back as NaN values.
    """
    padded_array = np.full(original_shape, np.nan)
    mask = np.ones(original_shape[0], dtype=bool)
    mask[removed_indices] = False
    padded_array[mask] = array
    return padded_array

def rotate_back(array, angle):
    """
    Rotates the array back to its original orientation.

    Parameters:
    - array (numpy.ndarray): The input array to be rotated back.
    - angle (float): The rotation angle in degrees.

    Returns:
    - rotated_back_array (numpy.ndarray): The array rotated back to its original orientation.
    """
    rotated_back_array = rotate(array, -angle, reshape=False, order=1, mode='constant', cval=np.nan)
    return rotated_back_array

A script and input file to test the functionality for rotating and rotating back can be found in
aeolis_rotation proposal.zip