Module geosardine.interpolate

@singledispatch
def idw(
    points: Union[str, np.ndarray],
    value: np.ndarray,
    spatial_res: Tuple[float, float],
    epsg: int = 4326,
    column_name: Optional[str] = None,
    longlat_distance: str = "harvesine",
    extent: Optional[Tuple[float, float, float, float]] = None,
    power: Union[float, int] = 2,
    distance_limit: float = 0.0,
) -> Optional[Raster]:
    """
    create interpolated raster from point by using Inverse Distance Weighting (Shepard)

    Parameters
    ----------
    points : numpy array, str
        list of points coordinate as numpy array or address of vector file  <br/>
        i.e shapefile or geojson  <br/>
        * if numpy array, then value input needed  <br/>
        * if str, then value is not needed instead will be created from file
    value : numpy array
        list of points value as numpy array, not needed if vector file used as input
    spatial_res : tuple or list of float
        spatial resolution in x and y axis
    column_name : str, default None
        column name needed to obtain value from attribute data of vector file  <br/>
        * If str, value will be read from respective column name  <br/>
        * If None, first column will be used as value
    epsg : int, default 4326
        EPSG code of reference system  <br/>
        * If 4326, WGS 1984 geographic system  <br/>
        * If int, epsg will be parsed
    longlat_distance: str harvesine or vincenty, default harvesine
        method used to calculate distance in spherical / ellipsoidal  <br/>
        * If harvesine, calculation will be faster but less accurate  <br/>
        * If vincenty, calculation will be slower but more accurate
    extent: tuple of float, default None
        how wide the raster will be  <br/>
        * If None, extent will be calculated from points input  <br/>
        * If tuple of float, user input of extent will be used
    power: float, default 2
        how smooth the interpolation will be
    distance_limit: float, default 0
        maximum distance to be interpolated, can't be negative

    Returns
    -------
    InterpolationResult

    Examples
    --------

    >>> xy = np.array([[106.8358,  -6.585 ],
    ...     [106.6039,  -6.7226],
    ...     [106.7589,  -6.4053],
    ...     [106.9674,  -6.7092],
    ...     [106.7956,  -6.5988]
    ... ])

    >>> values = np.array([132., 127.,  37.,  90., 182.])

    >>> idw(xy, values, spatial_res=(0.01,0.01), epsg=4326)

    >>> print(interpolated.array)
    [[ 88.63769859  86.24219616  83.60463194 ... 101.98185127 103.37001289 104.54621272]
     [ 90.12053232  87.79279317  85.22030848 ... 103.77118852 105.01425289 106.05302554]
     [ 91.82987695  89.60855271  87.14722258 ... 105.70090081 106.76928067 107.64635337]
     ...
     [127.21214817 127.33208302 127.53878268 ...  97.80436475  94.96247196 93.12113458]
     [127.11315081 127.18465002 127.33444124 ...  95.86455668  93.19212577 91.51135399]
     [127.0435062  127.0827023  127.19214624 ...  94.80175756  92.30685734 90.75707134]]
    """
    print("only support numpy array or vector file such as shapefile and geojson")
    return None

def idw_single(
    point: List[float],
    known_coordinates: np.ndarray,
    known_value: np.ndarray,
    epsg: int = 4326,
    longlat_distance: str = "harvesine",
    power: Union[float, int] = 2,
    distance_limit: float = 0.0,
) -> float:
    """

    Parameters
    ----------
    point : list
        list of single point to be interpolated
    known_coordinates : numpy array
        list of points coordinate as numpy array
    known_value: numpy array
        list of points value as numpy array, not needed if vector file used as input
    epsg : int, default 4326
        EPSG code of reference system  <br/>
        * If 4326, WGS 1984 geographic system  <br/>
        * If int, epsg will be parsed
    longlat_distance: str harvesine or vincenty, default harvesine
        method used to calculate distance in spherical / ellipsoidal  <br/>
        * If harvesine, calculation will be faster but less accurate  <br/>
        * If vincenty, calculation will be slower but more accurate
    power: float, default 2
        how smooth the interpolation will be
    distance_limit: float, default 0
        maximum distance to be interpolated, can't be negative

    Returns
    -------
    float
        interpolated value

    Examples
    --------
    >>> from geosardine.interpolate import idw_single

    >>> result = idw_single(
    ...     [860209, 9295740],
    ...     np.array([[767984, 9261620], [838926, 9234594]]),
    ...     np.array([[101.1, 102.2]]),
    ...     epsg=32748,
    ...     distance_limit=0
    ... )

    >>> print(result)
    101.86735169471324

    """
    if len(point) > 2:
        raise ValueError("only for single point, input can't be more than 2 items")
    crs = CRS.from_epsg(epsg)
    distance_calculation = longlat_distance
    if crs.is_projected:
        distance_calculation = "projected"

    interpolated = _idw(
        known_coordinates,
        known_value,
        np.array([point]),
        calc_distance[distance_calculation],
        power,
        distance_limit,
    )
    return interpolated[0]

class InterpolationResult:
    """
    Class to interpolation result

    Attributes
    ----------
    array : numpy array
        array of interpolated value.
    coordinates : numpy array
        coordinate array of interpolated value  <br/>
        each pixel / grid is x and y or longitude and latitude  <br/>
    crs : `rasterio.crs.CRS`
        crs of interpolated value
    extent : tuple
        extent of interpolated  <br/>
        x min, y min, x max, y max
    source : pathlib.Path, default None
        source file location  <br/>
        * if None, there is no source file  <br/>
        * if str, location of point file
    """

    def __init__(
        self,
        array: np.ndarray,
        coordinates: np.ndarray,
        crs: CRS,
        extent: Optional[Tuple[float, float, float, float]] = None,
        source: Optional[Union[str, Path]] = None,
    ):
        """
        Constructs interpolation result
        Parameters
        ----------
        array : numpy array
            array of interpolated value
        coordinates : numpy array
            coordinate array of interpolated value <br/>
            each pixel / grid is x and y or longitude and latitude
        crs : `rasterio.crs.CRS`
            crs of interpolated value
        extent : tuple, default None
            extent of interpolated <br/>
            x min, y min, x max, y max <br/>
            * if None, extent will be calculated from coordinate <br/>
            * if tuple, extent will be same as input
        source : str, pathlib.Path, default None
            source file location <br/>
            * if None, there is no source file <br/>
            * if str or `pathlib.Path`, location of point file
        """
        self.array = array
        self.crs = crs
        self.extent = extent
        self.source = source
        self.output: Optional[Path] = None
        self.affine: Affine = calc_affine(coordinates)
        if extent is None:
            self.extent = calc_extent(coordinates)
        del coordinates

        if type(source) == str:
            self.source = Path(source)

    def save(self, location: Optional[Union[str, Path]] = None) -> None:
        """
        save interpolated array as geotif
        Parameters
        ----------
        location : str, pathlib.Path, default None
            output location  <br/>
            * if None, tiff will be saved in the same directory and same name as source
                will only work if source is not None  <br/>
            * if str or pathlib.Path, tiff will be saved in there  <br/>

        """
        if self.source is None and location is None:
            raise ValueError("Please provide output location")

        self.output = location
        if self.source is not None and location is None:
            self.output = self.source.parent / f"{self.source.stem}.tif"
            print(f"OUTPUT is not specified, will be saved as {self.output}")

        save_raster(self.output, self.array, affine=self.affine, crs=self.crs)