Incorrect topographical data retrievation

[pragayshourya]

Bug Description

I believe the current method for retrieving topographical features is incorrect because it is not accounting for the projection in the OGGM module.

In the existing code:

def _retrieve_topo_features(
    df: pd.DataFrame,
    glacier_directories: list,
    gdirs_gridded: list,
    grouped_stakes: pd.api.typing.DataFrameGroupBy,
    voi: list,
) -> None:
    """Find the nearest recorded point with topographical features on the glacier for each stake."""

    for gdir, gdir_grid in zip(glacier_directories, gdirs_gridded):
        lat = grouped_stakes.get_group(gdir.rgi_id)[["POINT_LAT"]].values.flatten()
        lon = grouped_stakes.get_group(gdir.rgi_id)[["POINT_LON"]].values.flatten()

        topo_data = (gdir_grid.sel(
            x=lon, y=lat,
            method="nearest")[voi].to_dataframe().reset_index(drop=True))

        df.loc[df["RGIId"] == gdir.rgi_id, voi] = topo_data[voi]

where, the gdirs_gridded refers to the respective glacier data. OGGM projects the topographical data into a new projection specified by a proj statement that can be read by projection libraries. Consequently, some values are NaN even though the points fall within the glaciers.

This issue can be addressed by transforming the points into the new coordinate system using GeoPandas as follows:

def _process_gdir(gdir, gdir_gridded):
    temp_gdf = self.df_coords[self.df_coords['RGIId'] == gdir.rgi_id].copy()
    temp_gdf.to_crs(proj.CRS.from_proj4(gdir_gridded.attrs['pyproj_srs']), inplace=True)

    x = xr.DataArray(temp_gdf.geometry.x, dims='location')
    y = xr.DataArray(temp_gdf.geometry.y, dims='location')
    gdir_gridded['new_ID'] = (['location'], temp_gdf.index.values)

    topo_data_list = (
        gdir_gridded.sel(x=x, y=y, method='nearest')[voi + ['new_ID']].to_dataframe().reset_index(drop=True)
    )
    return topo_data_list

[pragayshourya]

Also, when extracting climate variables, following code has been utilised (an example):

x = xr.dataarray(temp_gdf.geometry.x, dims="location")

but for topographical feature extraction following has been used:

topo_data = (gdir_grid.sel(
            x=lon, y=lat,
            method="nearest")[voi].to_dataframe().reset_index(drop=True))

xarray.sel(x = [a list of xs] and y = [a list of ys]) assumes that these are the x and y coordinates of a grid x by y, rather than point coordinates. Since we are assigning values based on point coordinates, it would make more sense (In my opinion) to extract points using DataArrays.