Warn.zeropadding for small islands

[ThomasRoosli]

Describe the bug
The function zeropadding "Produces a rectangular shaped map from a non-rectangular map (e.g., country). For this,
a regular grid is created in the rectangle enclosing the non-rectangular map." This fails for very small islands like Bermuda.

To Reproduce
Steps to reproduce the behavior/error:

1. take the points of the LitPop Bermuda exposure and try to create a rectangular grid using zeropadding
2. It creates a ValueError
3. It works for most/all other LitPop Exposures (e.g. Liechtenstein)

Code example:

from climada_petals.engine.warn import Warn
from climada.util.api_client import Client
client = Client()
exp = client.get_litpop('Bermuda')
# exp = client.get_litpop('Liechtenstein')
grid, coord_haiti = Warn.zeropadding(exp.gdf.latitude, exp.gdf.longitude, exp.gdf.value)

A possible fix is to do integer math based on arcmiliseconds (a solution developed at MeteoSwiss):
this code would need to be added to warn.py and to the Warn class definition.

from typing import List, Tuple, TypeVar
import numpy.typing as npt
T = TypeVar('T', bound=np.number)

ARC_MILLISECONDS_PER_DEGREE = 3_600_000

class Warn:
    @staticmethod
    def zeropadding(lat: npt.NDArray[T], lon: npt.NDArray[T], val: npt.NDArray[np.floating]):
        """Produces a rectangular shaped map from a non-rectangular map (e.g., country). For this,
        a regular grid is created in the rectangle enclosing the non-rectangular map. The values are
        mapped onto the regular grid according to their coordinates. The regular gird is filled with
        zeros where no values are defined.
        are defined. This only works if the lat lon values of the non-rectangular map can be
        accurately represented on a grid with a regular resolution.

        Parameters
        ----------
        lat : np.ndarray
            Latitudes of values of map.
        lon : np.ndarray
            Longitudes of values of map.
        val : np.ndarray
            Values of quantity of interest at every coordinate given.

        Returns
        ----------
        map_rec : np.ndarray
            Rectangular map with a value for every grid point. Padded with zeros where no values in
            input map.
        coord_rec : np.ndarray
            Longitudes and Latitudes of every value of the map.
        """
        if not (len(lat) == len(lon) == len(val)):
            raise AssertionError(f'Lengths of lat, lon and val are not equal: {len(lat)}, {len(lon)}, {len(val)}')

        # enforce using integer values during zero-padding
        # convert latitude and longitude values to arc milliseconds if no integer values are passed
        floating_coordinates = np.issubdtype(lat.dtype, np.floating) or np.issubdtype(lon.dtype, np.floating)
        if floating_coordinates:
            lat = to_arc_milliseconds(lat)
            lon = to_arc_milliseconds(lon)

        lat_min = min(lat)
        lat_max = max(lat)
        lon_min = min(lon)
        lon_max = max(lon)

        grid_res_lon = Warn._get_resolution(lon)
        grid_res_lat = Warn._get_resolution(lat)

        un_lat = np.arange(
            lat_min,
            lat_max+grid_res_lat,
            grid_res_lat
        )
        un_lon = np.arange(
            lon_min,
            lon_max+grid_res_lon,
            grid_res_lon
        )
        # convert back to floating point number if originally passed
        if floating_coordinates:
            un_lat = from_arc_milliseconds(un_lat, decimals=12)
            un_lon = from_arc_milliseconds(un_lon, decimals=12)

        i = ((lat - lat_min) / grid_res_lat).astype(int)
        j = ((lon - lon_min) / grid_res_lon).astype(int)
        map_rec = np.zeros((len(un_lat), len(un_lon)))
        map_rec[i, j] = val

        lons, lats = np.meshgrid(
            un_lon, un_lat
        )
        coord_rec = np.vstack((lats.flatten(), lons.flatten())).transpose()

        return map_rec, coord_rec

    @staticmethod
    def _get_resolution(values: np.ndarray[int]) -> int:
        if len(values) <= 1 or values.min() == values.max():
            return 1
        diffs = abs(np.diff(np.unique(values)))
        if len(diffs) == 1:
            return diffs[0]
        return np.gcd.reduce(diffs)


def to_arc_milliseconds(values: npt.NDArray[np.floating]) -> npt.NDArray[np.signedinteger]:
    return np.round(np.multiply(values, ARC_MILLISECONDS_PER_DEGREE)).astype(int)


def from_arc_milliseconds(values: npt.NDArray[np.signedinteger], decimals: int = 12)-> npt.NDArray[np.floating]:
    return np.round(np.divide(values, ARC_MILLISECONDS_PER_DEGREE), decimals=decimals)