From edff40d3b544ee06e82cda41d1dc22c9dcd6ba73 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Fri, 27 Sep 2024 15:38:06 +0100 Subject: [PATCH 1/9] copy files over ocf_datapipes --- ocf_data_sampler/select/select_spatial_slice.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/ocf_data_sampler/select/select_spatial_slice.py b/ocf_data_sampler/select/select_spatial_slice.py index 0e2f47e..b95ecec 100644 --- a/ocf_data_sampler/select/select_spatial_slice.py +++ b/ocf_data_sampler/select/select_spatial_slice.py @@ -5,15 +5,15 @@ import numpy as np import xarray as xr -from ocf_datapipes.utils import Location -from ocf_datapipes.utils.geospatial import ( +from ocf_data_sampler.utils.location import Location +from ocf_data_sampler.utils.geospatial import ( lon_lat_to_geostationary_area_coords, lon_lat_to_osgb, osgb_to_geostationary_area_coords, osgb_to_lon_lat, spatial_coord_type, ) -from ocf_datapipes.utils.utils import searchsorted + logger = logging.getLogger(__name__) @@ -130,13 +130,8 @@ def _get_idx_of_pixel_closest_to_poi_geostationary( f"{y} is not in the interval {da[y_dim].min().values}: {da[y_dim].max().values}" # Get the index into x and y nearest to x_center_geostationary and y_center_geostationary: - x_index_at_center = searchsorted( - da[x_dim].values, center_geostationary.x, assume_ascending=True - ) - - y_index_at_center = searchsorted( - da[y_dim].values, center_geostationary.y, assume_ascending=True - ) + x_index_at_center = np.searchsorted(da[x_dim].values, center_geostationary.x) + y_index_at_center = np.searchsorted(da[y_dim].values, center_geostationary.y) return Location(x=x_index_at_center, y=y_index_at_center, coordinate_system="idx") From 46858f194f7b7a58929251587753f1576c9d6b05 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Fri, 27 Sep 2024 15:50:31 +0100 Subject: [PATCH 2/9] add files --- ocf_data_sampler/utils/__init__.py | 0 ocf_data_sampler/utils/geospatial.py | 283 +++++++++++++++++++++++++++ ocf_data_sampler/utils/location.py | 63 ++++++ 3 files changed, 346 insertions(+) create mode 100644 ocf_data_sampler/utils/__init__.py create mode 100644 ocf_data_sampler/utils/geospatial.py create mode 100644 ocf_data_sampler/utils/location.py diff --git a/ocf_data_sampler/utils/__init__.py b/ocf_data_sampler/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ocf_data_sampler/utils/geospatial.py b/ocf_data_sampler/utils/geospatial.py new file mode 100644 index 0000000..ce57ab6 --- /dev/null +++ b/ocf_data_sampler/utils/geospatial.py @@ -0,0 +1,283 @@ +"""Geospatial functions""" + +from datetime import datetime +from numbers import Number +from typing import Union + +import numpy as np +import pandas as pd +import pvlib +import pyproj +import xarray as xr + +# OSGB is also called "OSGB 1936 / British National Grid -- United +# Kingdom Ordnance Survey". OSGB is used in many UK electricity +# system maps, and is used by the UK Met Office UKV model. OSGB is a +# Transverse Mercator projection, using 'easting' and 'northing' +# coordinates which are in meters. See https://epsg.io/27700 +OSGB36 = 27700 + +# WGS84 is short for "World Geodetic System 1984", used in GPS. Uses +# latitude and longitude. +WGS84 = 4326 + + +_osgb_to_lon_lat = pyproj.Transformer.from_crs( + crs_from=OSGB36, crs_to=WGS84, always_xy=True +).transform +_lon_lat_to_osgb = pyproj.Transformer.from_crs( + crs_from=WGS84, crs_to=OSGB36, always_xy=True +).transform +_geod = pyproj.Geod(ellps="WGS84") + + +def osgb_to_lon_lat( + x: Union[Number, np.ndarray], y: Union[Number, np.ndarray] +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Change OSGB coordinates to lon, lat. + + Args: + x: osgb east-west + y: osgb north-south + Return: 2-tuple of longitude (east-west), latitude (north-south) + """ + return _osgb_to_lon_lat(xx=x, yy=y) + + +def lon_lat_to_osgb( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Change lon-lat coordinates to OSGB. + + Args: + x: longitude east-west + y: latitude north-south + + Return: 2-tuple of OSGB x, y + """ + return _lon_lat_to_osgb(xx=x, yy=y) + + +def lon_lat_to_geostationary_area_coords( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], + xr_data: Union[xr.Dataset, xr.DataArray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Loads geostationary area and change from lon-lat to geostationaery coords + + Args: + x: Longitude east-west + y: Latitude north-south + xr_data: xarray object with geostationary area + + Returns: + Geostationary coords: x, y + """ + # Only load these if using geostationary projection + import pyproj + import pyresample + + try: + area_definition_yaml = xr_data.attrs["area"] + except KeyError: + area_definition_yaml = xr_data.data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( + area_definition_yaml + ) + geostationary_crs = geostationary_area_definition.crs + lonlat_to_geostationary = pyproj.Transformer.from_crs( + crs_from=WGS84, + crs_to=geostationary_crs, + always_xy=True, + ).transform + return lonlat_to_geostationary(xx=x, yy=y) + + +def osgb_to_geostationary_area_coords( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], + xr_data: Union[xr.Dataset, xr.DataArray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Loads geostationary area and transformation from OSGB to geostationary coords + + Args: + x: osgb east-west + y: osgb north-south + xr_data: xarray object with geostationary area + + Returns: + Geostationary coords: x, y + """ + # Only load these if using geostationary projection + import pyproj + import pyresample + + try: + area_definition_yaml = xr_data.attrs["area"] + except KeyError: + area_definition_yaml = xr_data.data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( + area_definition_yaml + ) + geostationary_crs = geostationary_area_definition.crs + osgb_to_geostationary = pyproj.Transformer.from_crs( + crs_from=OSGB36, crs_to=geostationary_crs, always_xy=True + ).transform + return osgb_to_geostationary(xx=x, yy=y) + + +def geostationary_area_coords_to_osgb( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], + xr_data: Union[xr.Dataset, xr.DataArray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Loads geostationary area and change from geostationary coords to OSGB + + Args: + x: geostationary x coord + y: geostationary y coord + xr_data: xarray object with geostationary area + + Returns: + OSGB x, OSGB y + """ + # Only load these if using geostationary projection + import pyproj + import pyresample + + try: + area_definition_yaml = xr_data.attrs["area"] + except KeyError: + area_definition_yaml = xr_data.data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( + area_definition_yaml + ) + geostationary_crs = geostationary_area_definition.crs + geostationary_to_osgb = pyproj.Transformer.from_crs( + crs_from=geostationary_crs, crs_to=OSGB36, always_xy=True + ).transform + return geostationary_to_osgb(xx=x, yy=y) + + +def geostationary_area_coords_to_lonlat( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], + xr_data: Union[xr.Dataset, xr.DataArray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Loads geostationary area and change from geostationary to lon-lat coords + + Args: + x: geostationary x coord + y: geostationary y coord + xr_data: xarray object with geostationary area + + Returns: + longitude, latitude + """ + # Only load these if using geostationary projection + import pyproj + import pyresample + + try: + area_definition_yaml = xr_data.attrs["area"] + except KeyError: + area_definition_yaml = xr_data.data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( + area_definition_yaml + ) + geostationary_crs = geostationary_area_definition.crs + geostationary_to_lonlat = pyproj.Transformer.from_crs( + crs_from=geostationary_crs, crs_to=WGS84, always_xy=True + ).transform + return geostationary_to_lonlat(xx=x, yy=y) + + +def calculate_azimuth_and_elevation_angle( + latitude: float, longitude: float, datestamps: list[datetime] +) -> pd.DataFrame: + """ + Calculation the azimuth angle, and the elevation angle for several datetamps. + + But for one specific lat/lon location + + More details see: + https://www.celestis.com/resources/faq/what-are-the-azimuth-and-elevation-of-a-satellite/ + + Args: + latitude: latitude of the pv site + longitude: longitude of the pv site + datestamps: list of datestamps to calculate the sun angles. i.e the sun moves from east to + west in the day. + + Returns: Pandas data frame with the index the same as 'datestamps', with columns of + "elevation" and "azimuth" that have been calculate. + + """ + # get the solor position + solpos = pvlib.solarposition.get_solarposition(datestamps, latitude, longitude) + + # extract the information we want + return solpos[["elevation", "azimuth"]] + + +def move_lon_lat_by_meters(lon, lat, meters_east, meters_north): + """ + Move a (lon, lat) by a certain number of meters north and east + + Args: + lon: longitude + lat: latitude + meters_east: number of meters to move east + meters_north: number of meters to move north + + Returns: + tuple of lon, lat + """ + new_lon = _geod.fwd(lons=lon, lats=lat, az=90, dist=meters_east)[0] + new_lat = _geod.fwd(lons=lon, lats=lat, az=0, dist=meters_north)[1] + return new_lon, new_lat + + +def _coord_priority(available_coords): + if "longitude" in available_coords: + return "lon_lat", "longitude", "latitude" + elif "x_geostationary" in available_coords: + return "geostationary", "x_geostationary", "y_geostationary" + elif "x_osgb" in available_coords: + return "osgb", "x_osgb", "y_osgb" + elif "x" in available_coords: + return "xy", "x", "y" + else: + return None, None, None + + +def spatial_coord_type(ds: xr.Dataset): + """Searches the dataset to determine the kind of spatial coordinates present. + + This search has a preference for the dimension coordinates of the xarray object. If none of the + expected coordinates exist in the dimension coordinates, it then searches the non-dimension + coordinates. See https://docs.xarray.dev/en/latest/user-guide/data-structures.html#coordinates. + + Args: + ds: Dataset with spatial coords + + Returns: + str: The kind of the coordinate system + x_coord: Name of the x-coordinate + y_coord: Name of the y-coordinate + """ + if isinstance(ds, xr.DataArray): + # Search dimension coords of dataarray + coords = _coord_priority(ds.xindexes) + elif isinstance(ds, xr.Dataset): + # Search dimension coords of all variables in dataset + coords = _coord_priority(set([v for k in ds.keys() for v in list(ds[k].xindexes)])) + else: + raise ValueError(f"Unrecognized input type: {type(ds)}") + + if coords == (None, None, None): + # If no dimension coords found, search non-dimension coords + coords = _coord_priority(list(ds.coords)) + + return coords diff --git a/ocf_data_sampler/utils/location.py b/ocf_data_sampler/utils/location.py new file mode 100644 index 0000000..fefb332 --- /dev/null +++ b/ocf_data_sampler/utils/location.py @@ -0,0 +1,63 @@ +"""location""" + +from typing import Optional + +import numpy as np +from pydantic import BaseModel, Field, validator + + +class Location(BaseModel): + """Represent a spatial location.""" + + coordinate_system: Optional[str] = "osgb" # ["osgb", "lon_lat", "geostationary", "idx"] + x: float + y: float + id: Optional[int] = Field(None) + + @validator("coordinate_system", pre=True, always=True) + def validate_coordinate_system(cls, v): + """Validate 'coordinate_system'""" + allowed_coordinate_systen = ["osgb", "lon_lat", "geostationary", "idx"] + if v not in allowed_coordinate_systen: + raise ValueError(f"coordinate_system = {v} is not in {allowed_coordinate_systen}") + return v + + @validator("x") + def validate_x(cls, v, values): + """Validate 'x'""" + min_x: float + max_x: float + if "coordinate_system" not in values: + raise ValueError("coordinate_system is incorrect") + co = values["coordinate_system"] + if co == "osgb": + min_x, max_x = -103976.3, 652897.98 + if co == "lon_lat": + min_x, max_x = -180, 180 + if co == "geostationary": + min_x, max_x = -5568748.275756836, 5567248.074173927 + if co == "idx": + min_x, max_x = 0, np.inf + if v < min_x or v > max_x: + raise ValueError(f"x = {v} must be within {[min_x, max_x]} for {co} coordinate system") + return v + + @validator("y") + def validate_y(cls, v, values): + """Validate 'y'""" + min_y: float + max_y: float + if "coordinate_system" not in values: + raise ValueError("coordinate_system is incorrect") + co = values["coordinate_system"] + if co == "osgb": + min_y, max_y = -16703.87, 1199851.44 + if co == "lon_lat": + min_y, max_y = -90, 90 + if co == "geostationary": + min_y, max_y = 1393687.2151494026, 5570748.323202133 + if co == "idx": + min_y, max_y = 0, np.inf + if v < min_y or v > max_y: + raise ValueError(f"y = {v} must be within {[min_y, max_y]} for {co} coordinate system") + return v From b42da3d3a8004099f4eda7cbedf032cb0aee8cdd Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Fri, 27 Sep 2024 16:17:37 +0100 Subject: [PATCH 3/9] refactor --- ocf_data_sampler/select/geospatial.py | 137 +++++++++ ocf_data_sampler/select/location.py | 62 ++++ .../select/select_spatial_slice.py | 12 +- ocf_data_sampler/utils/__init__.py | 0 ocf_data_sampler/utils/geospatial.py | 283 ------------------ ocf_data_sampler/utils/location.py | 63 ---- 6 files changed, 203 insertions(+), 354 deletions(-) create mode 100644 ocf_data_sampler/select/geospatial.py create mode 100644 ocf_data_sampler/select/location.py delete mode 100644 ocf_data_sampler/utils/__init__.py delete mode 100644 ocf_data_sampler/utils/geospatial.py delete mode 100644 ocf_data_sampler/utils/location.py diff --git a/ocf_data_sampler/select/geospatial.py b/ocf_data_sampler/select/geospatial.py new file mode 100644 index 0000000..6be26b3 --- /dev/null +++ b/ocf_data_sampler/select/geospatial.py @@ -0,0 +1,137 @@ +"""Geospatial functions""" + +from datetime import datetime +from numbers import Number +from typing import Union + +import numpy as np +import pandas as pd +import pvlib +import pyproj +import xarray as xr + +# OSGB is also called "OSGB 1936 / British National Grid -- United +# Kingdom Ordnance Survey". OSGB is used in many UK electricity +# system maps, and is used by the UK Met Office UKV model. OSGB is a +# Transverse Mercator projection, using 'easting' and 'northing' +# coordinates which are in meters. See https://epsg.io/27700 +OSGB36 = 27700 + +# WGS84 is short for "World Geodetic System 1984", used in GPS. Uses +# latitude and longitude. +WGS84 = 4326 + + +_osgb_to_lon_lat = pyproj.Transformer.from_crs( + crs_from=OSGB36, crs_to=WGS84, always_xy=True +).transform +_lon_lat_to_osgb = pyproj.Transformer.from_crs( + crs_from=WGS84, crs_to=OSGB36, always_xy=True +).transform +_geod = pyproj.Geod(ellps="WGS84") + + +def osgb_to_lon_lat( + x: Union[Number, np.ndarray], y: Union[Number, np.ndarray] +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Change OSGB coordinates to lon, lat. + + Args: + x: osgb east-west + y: osgb north-south + Return: 2-tuple of longitude (east-west), latitude (north-south) + """ + return _osgb_to_lon_lat(xx=x, yy=y) + + +def lon_lat_to_osgb( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Change lon-lat coordinates to OSGB. + + Args: + x: longitude east-west + y: latitude north-south + + Return: 2-tuple of OSGB x, y + """ + return _lon_lat_to_osgb(xx=x, yy=y) + + + +def osgb_to_geostationary_area_coords( + x: Union[Number, np.ndarray], + y: Union[Number, np.ndarray], + xr_data: Union[xr.Dataset, xr.DataArray], +) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: + """Loads geostationary area and transformation from OSGB to geostationary coords + + Args: + x: osgb east-west + y: osgb north-south + xr_data: xarray object with geostationary area + + Returns: + Geostationary coords: x, y + """ + # Only load these if using geostationary projection + import pyproj + import pyresample + + try: + area_definition_yaml = xr_data.attrs["area"] + except KeyError: + area_definition_yaml = xr_data.data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( + area_definition_yaml + ) + geostationary_crs = geostationary_area_definition.crs + osgb_to_geostationary = pyproj.Transformer.from_crs( + crs_from=OSGB36, crs_to=geostationary_crs, always_xy=True + ).transform + return osgb_to_geostationary(xx=x, yy=y) + + +def _coord_priority(available_coords): + if "longitude" in available_coords: + return "lon_lat", "longitude", "latitude" + elif "x_geostationary" in available_coords: + return "geostationary", "x_geostationary", "y_geostationary" + elif "x_osgb" in available_coords: + return "osgb", "x_osgb", "y_osgb" + elif "x" in available_coords: + return "xy", "x", "y" + else: + return None, None, None + + +def spatial_coord_type(ds: xr.Dataset): + """Searches the dataset to determine the kind of spatial coordinates present. + + This search has a preference for the dimension coordinates of the xarray object. If none of the + expected coordinates exist in the dimension coordinates, it then searches the non-dimension + coordinates. See https://docs.xarray.dev/en/latest/user-guide/data-structures.html#coordinates. + + Args: + ds: Dataset with spatial coords + + Returns: + str: The kind of the coordinate system + x_coord: Name of the x-coordinate + y_coord: Name of the y-coordinate + """ + if isinstance(ds, xr.DataArray): + # Search dimension coords of dataarray + coords = _coord_priority(ds.xindexes) + elif isinstance(ds, xr.Dataset): + # Search dimension coords of all variables in dataset + coords = _coord_priority(set([v for k in ds.keys() for v in list(ds[k].xindexes)])) + else: + raise ValueError(f"Unrecognized input type: {type(ds)}") + + if coords == (None, None, None): + # If no dimension coords found, search non-dimension coords + coords = _coord_priority(list(ds.coords)) + + return coords diff --git a/ocf_data_sampler/select/location.py b/ocf_data_sampler/select/location.py new file mode 100644 index 0000000..9cfa9cf --- /dev/null +++ b/ocf_data_sampler/select/location.py @@ -0,0 +1,62 @@ +"""location""" + +from typing import Optional + +import numpy as np +from pydantic import BaseModel, Field, model_validator + + +allowed_coordinate_systems =["osgb", "lon_lat", "geostationary", "idx"] + +class Location(BaseModel): + """Represent a spatial location.""" + + coordinate_system: Optional[str] = "osgb" # ["osgb", "lon_lat", "geostationary", "idx"] + x: float + y: float + id: Optional[int] = Field(None) + + @model_validator(mode='after') + def validate_coordinate_system(self): + """Validate 'coordinate_system'""" + if self.coordinate_system not in allowed_coordinate_systems: + raise ValueError(f"coordinate_system = {self.coordinate_system} is not in {allowed_coordinate_systems}") + return self + + @model_validator(mode='after') + def validate_x(self): + """Validate 'x'""" + min_x: float + max_x: float + + co = self.coordinate_system + if co == "osgb": + min_x, max_x = -103976.3, 652897.98 + if co == "lon_lat": + min_x, max_x = -180, 180 + if co == "geostationary": + min_x, max_x = -5568748.275756836, 5567248.074173927 + if co == "idx": + min_x, max_x = 0, np.inf + if self.x < min_x or self.x > max_x: + raise ValueError(f"x = {self.x} must be within {[min_x, max_x]} for {co} coordinate system") + return self + + @model_validator(mode='after') + def validate_y(self): + """Validate 'y'""" + min_y: float + max_y: float + + co = self.coordinate_system + if co == "osgb": + min_y, max_y = -16703.87, 1199851.44 + if co == "lon_lat": + min_y, max_y = -90, 90 + if co == "geostationary": + min_y, max_y = 1393687.2151494026, 5570748.323202133 + if co == "idx": + min_y, max_y = 0, np.inf + if self.y < min_y or self.y > max_y: + raise ValueError(f"y = {self.y} must be within {[min_y, max_y]} for {co} coordinate system") + return self diff --git a/ocf_data_sampler/select/select_spatial_slice.py b/ocf_data_sampler/select/select_spatial_slice.py index b95ecec..f1b2ead 100644 --- a/ocf_data_sampler/select/select_spatial_slice.py +++ b/ocf_data_sampler/select/select_spatial_slice.py @@ -5,9 +5,8 @@ import numpy as np import xarray as xr -from ocf_data_sampler.utils.location import Location -from ocf_data_sampler.utils.geospatial import ( - lon_lat_to_geostationary_area_coords, +from ocf_data_sampler.select.location import Location +from ocf_data_sampler.select.geospatial import ( lon_lat_to_osgb, osgb_to_geostationary_area_coords, osgb_to_lon_lat, @@ -45,9 +44,6 @@ def convert_coords_to_match_xarray( if from_coords == "osgb": x, y = osgb_to_geostationary_area_coords(x, y, da) - elif from_coords == "lon_lat": - x, y = lon_lat_to_geostationary_area_coords(x, y, da) - elif target_coords == "lon_lat": if from_coords == "osgb": x, y = osgb_to_lon_lat(x, y) @@ -97,8 +93,8 @@ def _get_idx_of_pixel_closest_to_poi( x_index = da.get_index(x_dim) y_index = da.get_index(y_dim) - closest_x = x_index.get_indexer([x], method="nearest")[0] - closest_y = y_index.get_indexer([y], method="nearest")[0] + closest_x = float(x_index.get_indexer([x], method="nearest")[0]) + closest_y = float(y_index.get_indexer([y], method="nearest")[0]) return Location(x=closest_x, y=closest_y, coordinate_system="idx") diff --git a/ocf_data_sampler/utils/__init__.py b/ocf_data_sampler/utils/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/ocf_data_sampler/utils/geospatial.py b/ocf_data_sampler/utils/geospatial.py deleted file mode 100644 index ce57ab6..0000000 --- a/ocf_data_sampler/utils/geospatial.py +++ /dev/null @@ -1,283 +0,0 @@ -"""Geospatial functions""" - -from datetime import datetime -from numbers import Number -from typing import Union - -import numpy as np -import pandas as pd -import pvlib -import pyproj -import xarray as xr - -# OSGB is also called "OSGB 1936 / British National Grid -- United -# Kingdom Ordnance Survey". OSGB is used in many UK electricity -# system maps, and is used by the UK Met Office UKV model. OSGB is a -# Transverse Mercator projection, using 'easting' and 'northing' -# coordinates which are in meters. See https://epsg.io/27700 -OSGB36 = 27700 - -# WGS84 is short for "World Geodetic System 1984", used in GPS. Uses -# latitude and longitude. -WGS84 = 4326 - - -_osgb_to_lon_lat = pyproj.Transformer.from_crs( - crs_from=OSGB36, crs_to=WGS84, always_xy=True -).transform -_lon_lat_to_osgb = pyproj.Transformer.from_crs( - crs_from=WGS84, crs_to=OSGB36, always_xy=True -).transform -_geod = pyproj.Geod(ellps="WGS84") - - -def osgb_to_lon_lat( - x: Union[Number, np.ndarray], y: Union[Number, np.ndarray] -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Change OSGB coordinates to lon, lat. - - Args: - x: osgb east-west - y: osgb north-south - Return: 2-tuple of longitude (east-west), latitude (north-south) - """ - return _osgb_to_lon_lat(xx=x, yy=y) - - -def lon_lat_to_osgb( - x: Union[Number, np.ndarray], - y: Union[Number, np.ndarray], -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Change lon-lat coordinates to OSGB. - - Args: - x: longitude east-west - y: latitude north-south - - Return: 2-tuple of OSGB x, y - """ - return _lon_lat_to_osgb(xx=x, yy=y) - - -def lon_lat_to_geostationary_area_coords( - x: Union[Number, np.ndarray], - y: Union[Number, np.ndarray], - xr_data: Union[xr.Dataset, xr.DataArray], -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Loads geostationary area and change from lon-lat to geostationaery coords - - Args: - x: Longitude east-west - y: Latitude north-south - xr_data: xarray object with geostationary area - - Returns: - Geostationary coords: x, y - """ - # Only load these if using geostationary projection - import pyproj - import pyresample - - try: - area_definition_yaml = xr_data.attrs["area"] - except KeyError: - area_definition_yaml = xr_data.data.attrs["area"] - geostationary_area_definition = pyresample.area_config.load_area_from_string( - area_definition_yaml - ) - geostationary_crs = geostationary_area_definition.crs - lonlat_to_geostationary = pyproj.Transformer.from_crs( - crs_from=WGS84, - crs_to=geostationary_crs, - always_xy=True, - ).transform - return lonlat_to_geostationary(xx=x, yy=y) - - -def osgb_to_geostationary_area_coords( - x: Union[Number, np.ndarray], - y: Union[Number, np.ndarray], - xr_data: Union[xr.Dataset, xr.DataArray], -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Loads geostationary area and transformation from OSGB to geostationary coords - - Args: - x: osgb east-west - y: osgb north-south - xr_data: xarray object with geostationary area - - Returns: - Geostationary coords: x, y - """ - # Only load these if using geostationary projection - import pyproj - import pyresample - - try: - area_definition_yaml = xr_data.attrs["area"] - except KeyError: - area_definition_yaml = xr_data.data.attrs["area"] - geostationary_area_definition = pyresample.area_config.load_area_from_string( - area_definition_yaml - ) - geostationary_crs = geostationary_area_definition.crs - osgb_to_geostationary = pyproj.Transformer.from_crs( - crs_from=OSGB36, crs_to=geostationary_crs, always_xy=True - ).transform - return osgb_to_geostationary(xx=x, yy=y) - - -def geostationary_area_coords_to_osgb( - x: Union[Number, np.ndarray], - y: Union[Number, np.ndarray], - xr_data: Union[xr.Dataset, xr.DataArray], -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Loads geostationary area and change from geostationary coords to OSGB - - Args: - x: geostationary x coord - y: geostationary y coord - xr_data: xarray object with geostationary area - - Returns: - OSGB x, OSGB y - """ - # Only load these if using geostationary projection - import pyproj - import pyresample - - try: - area_definition_yaml = xr_data.attrs["area"] - except KeyError: - area_definition_yaml = xr_data.data.attrs["area"] - geostationary_area_definition = pyresample.area_config.load_area_from_string( - area_definition_yaml - ) - geostationary_crs = geostationary_area_definition.crs - geostationary_to_osgb = pyproj.Transformer.from_crs( - crs_from=geostationary_crs, crs_to=OSGB36, always_xy=True - ).transform - return geostationary_to_osgb(xx=x, yy=y) - - -def geostationary_area_coords_to_lonlat( - x: Union[Number, np.ndarray], - y: Union[Number, np.ndarray], - xr_data: Union[xr.Dataset, xr.DataArray], -) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: - """Loads geostationary area and change from geostationary to lon-lat coords - - Args: - x: geostationary x coord - y: geostationary y coord - xr_data: xarray object with geostationary area - - Returns: - longitude, latitude - """ - # Only load these if using geostationary projection - import pyproj - import pyresample - - try: - area_definition_yaml = xr_data.attrs["area"] - except KeyError: - area_definition_yaml = xr_data.data.attrs["area"] - geostationary_area_definition = pyresample.area_config.load_area_from_string( - area_definition_yaml - ) - geostationary_crs = geostationary_area_definition.crs - geostationary_to_lonlat = pyproj.Transformer.from_crs( - crs_from=geostationary_crs, crs_to=WGS84, always_xy=True - ).transform - return geostationary_to_lonlat(xx=x, yy=y) - - -def calculate_azimuth_and_elevation_angle( - latitude: float, longitude: float, datestamps: list[datetime] -) -> pd.DataFrame: - """ - Calculation the azimuth angle, and the elevation angle for several datetamps. - - But for one specific lat/lon location - - More details see: - https://www.celestis.com/resources/faq/what-are-the-azimuth-and-elevation-of-a-satellite/ - - Args: - latitude: latitude of the pv site - longitude: longitude of the pv site - datestamps: list of datestamps to calculate the sun angles. i.e the sun moves from east to - west in the day. - - Returns: Pandas data frame with the index the same as 'datestamps', with columns of - "elevation" and "azimuth" that have been calculate. - - """ - # get the solor position - solpos = pvlib.solarposition.get_solarposition(datestamps, latitude, longitude) - - # extract the information we want - return solpos[["elevation", "azimuth"]] - - -def move_lon_lat_by_meters(lon, lat, meters_east, meters_north): - """ - Move a (lon, lat) by a certain number of meters north and east - - Args: - lon: longitude - lat: latitude - meters_east: number of meters to move east - meters_north: number of meters to move north - - Returns: - tuple of lon, lat - """ - new_lon = _geod.fwd(lons=lon, lats=lat, az=90, dist=meters_east)[0] - new_lat = _geod.fwd(lons=lon, lats=lat, az=0, dist=meters_north)[1] - return new_lon, new_lat - - -def _coord_priority(available_coords): - if "longitude" in available_coords: - return "lon_lat", "longitude", "latitude" - elif "x_geostationary" in available_coords: - return "geostationary", "x_geostationary", "y_geostationary" - elif "x_osgb" in available_coords: - return "osgb", "x_osgb", "y_osgb" - elif "x" in available_coords: - return "xy", "x", "y" - else: - return None, None, None - - -def spatial_coord_type(ds: xr.Dataset): - """Searches the dataset to determine the kind of spatial coordinates present. - - This search has a preference for the dimension coordinates of the xarray object. If none of the - expected coordinates exist in the dimension coordinates, it then searches the non-dimension - coordinates. See https://docs.xarray.dev/en/latest/user-guide/data-structures.html#coordinates. - - Args: - ds: Dataset with spatial coords - - Returns: - str: The kind of the coordinate system - x_coord: Name of the x-coordinate - y_coord: Name of the y-coordinate - """ - if isinstance(ds, xr.DataArray): - # Search dimension coords of dataarray - coords = _coord_priority(ds.xindexes) - elif isinstance(ds, xr.Dataset): - # Search dimension coords of all variables in dataset - coords = _coord_priority(set([v for k in ds.keys() for v in list(ds[k].xindexes)])) - else: - raise ValueError(f"Unrecognized input type: {type(ds)}") - - if coords == (None, None, None): - # If no dimension coords found, search non-dimension coords - coords = _coord_priority(list(ds.coords)) - - return coords diff --git a/ocf_data_sampler/utils/location.py b/ocf_data_sampler/utils/location.py deleted file mode 100644 index fefb332..0000000 --- a/ocf_data_sampler/utils/location.py +++ /dev/null @@ -1,63 +0,0 @@ -"""location""" - -from typing import Optional - -import numpy as np -from pydantic import BaseModel, Field, validator - - -class Location(BaseModel): - """Represent a spatial location.""" - - coordinate_system: Optional[str] = "osgb" # ["osgb", "lon_lat", "geostationary", "idx"] - x: float - y: float - id: Optional[int] = Field(None) - - @validator("coordinate_system", pre=True, always=True) - def validate_coordinate_system(cls, v): - """Validate 'coordinate_system'""" - allowed_coordinate_systen = ["osgb", "lon_lat", "geostationary", "idx"] - if v not in allowed_coordinate_systen: - raise ValueError(f"coordinate_system = {v} is not in {allowed_coordinate_systen}") - return v - - @validator("x") - def validate_x(cls, v, values): - """Validate 'x'""" - min_x: float - max_x: float - if "coordinate_system" not in values: - raise ValueError("coordinate_system is incorrect") - co = values["coordinate_system"] - if co == "osgb": - min_x, max_x = -103976.3, 652897.98 - if co == "lon_lat": - min_x, max_x = -180, 180 - if co == "geostationary": - min_x, max_x = -5568748.275756836, 5567248.074173927 - if co == "idx": - min_x, max_x = 0, np.inf - if v < min_x or v > max_x: - raise ValueError(f"x = {v} must be within {[min_x, max_x]} for {co} coordinate system") - return v - - @validator("y") - def validate_y(cls, v, values): - """Validate 'y'""" - min_y: float - max_y: float - if "coordinate_system" not in values: - raise ValueError("coordinate_system is incorrect") - co = values["coordinate_system"] - if co == "osgb": - min_y, max_y = -16703.87, 1199851.44 - if co == "lon_lat": - min_y, max_y = -90, 90 - if co == "geostationary": - min_y, max_y = 1393687.2151494026, 5570748.323202133 - if co == "idx": - min_y, max_y = 0, np.inf - if v < min_y or v > max_y: - raise ValueError(f"y = {v} must be within {[min_y, max_y]} for {co} coordinate system") - return v From dd297b0ae50c136e1cacec9c1ec370aaf3c07e44 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Fri, 27 Sep 2024 16:26:46 +0100 Subject: [PATCH 4/9] tidy up --- ocf_data_sampler/select/geospatial.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/ocf_data_sampler/select/geospatial.py b/ocf_data_sampler/select/geospatial.py index 6be26b3..34149f5 100644 --- a/ocf_data_sampler/select/geospatial.py +++ b/ocf_data_sampler/select/geospatial.py @@ -59,7 +59,6 @@ def lon_lat_to_osgb( return _lon_lat_to_osgb(xx=x, yy=y) - def osgb_to_geostationary_area_coords( x: Union[Number, np.ndarray], y: Union[Number, np.ndarray], @@ -100,18 +99,14 @@ def _coord_priority(available_coords): return "geostationary", "x_geostationary", "y_geostationary" elif "x_osgb" in available_coords: return "osgb", "x_osgb", "y_osgb" - elif "x" in available_coords: - return "xy", "x", "y" else: - return None, None, None + raise ValueError(f"Unrecognized coordinate system: {available_coords}") def spatial_coord_type(ds: xr.Dataset): """Searches the dataset to determine the kind of spatial coordinates present. - This search has a preference for the dimension coordinates of the xarray object. If none of the - expected coordinates exist in the dimension coordinates, it then searches the non-dimension - coordinates. See https://docs.xarray.dev/en/latest/user-guide/data-structures.html#coordinates. + This search has a preference for the dimension coordinates of the xarray object. Args: ds: Dataset with spatial coords @@ -124,14 +119,7 @@ def spatial_coord_type(ds: xr.Dataset): if isinstance(ds, xr.DataArray): # Search dimension coords of dataarray coords = _coord_priority(ds.xindexes) - elif isinstance(ds, xr.Dataset): - # Search dimension coords of all variables in dataset - coords = _coord_priority(set([v for k in ds.keys() for v in list(ds[k].xindexes)])) else: raise ValueError(f"Unrecognized input type: {type(ds)}") - if coords == (None, None, None): - # If no dimension coords found, search non-dimension coords - coords = _coord_priority(list(ds.coords)) - return coords From ef7c83ee6b0f110cff07100821034979376c6dae Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Fri, 27 Sep 2024 16:28:47 +0100 Subject: [PATCH 5/9] att test location --- tests/select/test_location.py | 67 +++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 tests/select/test_location.py diff --git a/tests/select/test_location.py b/tests/select/test_location.py new file mode 100644 index 0000000..763f372 --- /dev/null +++ b/tests/select/test_location.py @@ -0,0 +1,67 @@ +from ocf_data_sampler.select.location import Location +import pytest + + +def test_make_valid_location_object_with_default_coordinate_system(): + x, y = -1000.5, 50000 + location = Location(x=x, y=y) + assert location.x == x, "location.x value not set correctly" + assert location.y == y, "location.x value not set correctly" + assert ( + location.coordinate_system == "osgb" + ), "location.coordinate_system value not set correctly" + + +def test_make_valid_location_object_with_osgb_coordinate_system(): + x, y, coordinate_system = 1.2, 22.9, "osgb" + location = Location(x=x, y=y, coordinate_system=coordinate_system) + assert location.x == x, "location.x value not set correctly" + assert location.y == y, "location.x value not set correctly" + assert ( + location.coordinate_system == coordinate_system + ), "location.coordinate_system value not set correctly" + + +def test_make_valid_location_object_with_lon_lat_coordinate_system(): + x, y, coordinate_system = 1.2, 1.2, "lon_lat" + location = Location(x=x, y=y, coordinate_system=coordinate_system) + assert location.x == x, "location.x value not set correctly" + assert location.y == y, "location.x value not set correctly" + assert ( + location.coordinate_system == coordinate_system + ), "location.coordinate_system value not set correctly" + + +def test_make_invalid_location_object_with_invalid_osgb_x(): + x, y, coordinate_system = 10000000, 1.2, "osgb" + with pytest.raises(ValueError) as err: + _ = Location(x=x, y=y, coordinate_system=coordinate_system) + assert err.typename == "ValidationError" + + +def test_make_invalid_location_object_with_invalid_osgb_y(): + x, y, coordinate_system = 2.5, 10000000, "osgb" + with pytest.raises(ValueError) as err: + _ = Location(x=x, y=y, coordinate_system=coordinate_system) + assert err.typename == "ValidationError" + + +def test_make_invalid_location_object_with_invalid_lon_lat_x(): + x, y, coordinate_system = 200, 1.2, "lon_lat" + with pytest.raises(ValueError) as err: + _ = Location(x=x, y=y, coordinate_system=coordinate_system) + assert err.typename == "ValidationError" + + +def test_make_invalid_location_object_with_invalid_lon_lat_y(): + x, y, coordinate_system = 2.5, -200, "lon_lat" + with pytest.raises(ValueError) as err: + _ = Location(x=x, y=y, coordinate_system=coordinate_system) + assert err.typename == "ValidationError" + + +def test_make_invalid_location_object_with_invalid_coordinate_system(): + x, y, coordinate_system = 2.5, 1000, "abcd" + with pytest.raises(ValueError) as err: + _ = Location(x=x, y=y, coordinate_system=coordinate_system) + assert err.typename == "ValidationError" From 92056a979198437132e7d2b239082dc1da54bd94 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Mon, 30 Sep 2024 14:10:37 +0100 Subject: [PATCH 6/9] tidy up imports, typing, PR comments --- ocf_data_sampler/select/geospatial.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/ocf_data_sampler/select/geospatial.py b/ocf_data_sampler/select/geospatial.py index 34149f5..fb15a8a 100644 --- a/ocf_data_sampler/select/geospatial.py +++ b/ocf_data_sampler/select/geospatial.py @@ -1,12 +1,9 @@ """Geospatial functions""" -from datetime import datetime from numbers import Number from typing import Union import numpy as np -import pandas as pd -import pvlib import pyproj import xarray as xr @@ -75,7 +72,6 @@ def osgb_to_geostationary_area_coords( Geostationary coords: x, y """ # Only load these if using geostationary projection - import pyproj import pyresample try: @@ -103,8 +99,8 @@ def _coord_priority(available_coords): raise ValueError(f"Unrecognized coordinate system: {available_coords}") -def spatial_coord_type(ds: xr.Dataset): - """Searches the dataset to determine the kind of spatial coordinates present. +def spatial_coord_type(ds: xr.DataArray): + """Searches the data array to determine the kind of spatial coordinates present. This search has a preference for the dimension coordinates of the xarray object. From cb2a8aa364a6954357e68068cdbd766ffff29fd5 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Mon, 30 Sep 2024 17:00:25 +0100 Subject: [PATCH 7/9] trim down for PR comments --- ocf_data_sampler/select/geospatial.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/ocf_data_sampler/select/geospatial.py b/ocf_data_sampler/select/geospatial.py index fb15a8a..3841395 100644 --- a/ocf_data_sampler/select/geospatial.py +++ b/ocf_data_sampler/select/geospatial.py @@ -59,7 +59,7 @@ def lon_lat_to_osgb( def osgb_to_geostationary_area_coords( x: Union[Number, np.ndarray], y: Union[Number, np.ndarray], - xr_data: Union[xr.Dataset, xr.DataArray], + xr_data: xr.DataArray, ) -> tuple[Union[Number, np.ndarray], Union[Number, np.ndarray]]: """Loads geostationary area and transformation from OSGB to geostationary coords @@ -74,10 +74,8 @@ def osgb_to_geostationary_area_coords( # Only load these if using geostationary projection import pyresample - try: - area_definition_yaml = xr_data.attrs["area"] - except KeyError: - area_definition_yaml = xr_data.data.attrs["area"] + area_definition_yaml = xr_data.attrs["area"] + geostationary_area_definition = pyresample.area_config.load_area_from_string( area_definition_yaml ) From 32697d8635ce5ee2df3f94f72e80152112c76931 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Mon, 30 Sep 2024 19:33:03 +0100 Subject: [PATCH 8/9] tidy up --- ocf_data_sampler/select/select_spatial_slice.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ocf_data_sampler/select/select_spatial_slice.py b/ocf_data_sampler/select/select_spatial_slice.py index f1b2ead..9dad8ea 100644 --- a/ocf_data_sampler/select/select_spatial_slice.py +++ b/ocf_data_sampler/select/select_spatial_slice.py @@ -93,8 +93,8 @@ def _get_idx_of_pixel_closest_to_poi( x_index = da.get_index(x_dim) y_index = da.get_index(y_dim) - closest_x = float(x_index.get_indexer([x], method="nearest")[0]) - closest_y = float(y_index.get_indexer([y], method="nearest")[0]) + closest_x = x_index.get_indexer([x], method="nearest")[0] + closest_y = y_index.get_indexer([y], method="nearest")[0] return Location(x=closest_x, y=closest_y, coordinate_system="idx") From 24d8162194734a2a05786b066176596d1c4f4c18 Mon Sep 17 00:00:00 2001 From: peterdudfield Date: Tue, 1 Oct 2024 11:18:37 +0100 Subject: [PATCH 9/9] strip out un needed --- ocf_data_sampler/select/geospatial.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ocf_data_sampler/select/geospatial.py b/ocf_data_sampler/select/geospatial.py index 3841395..8137f16 100644 --- a/ocf_data_sampler/select/geospatial.py +++ b/ocf_data_sampler/select/geospatial.py @@ -25,7 +25,6 @@ _lon_lat_to_osgb = pyproj.Transformer.from_crs( crs_from=WGS84, crs_to=OSGB36, always_xy=True ).transform -_geod = pyproj.Geod(ellps="WGS84") def osgb_to_lon_lat(