Feat/radiance reflectance and l457 (#9)

* Query landsat imagery and add function to obtain integrated bands from hyperspectral images

* Bug in figure out Landsat-7

* Added to download ts of landsat 4, 5 and 7

---------

Co-authored-by: Gonzalo Mateo Garcia <[email protected]>

georeader/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.1.3"
+__version__ = "1.1.4"
 
 import math
 from typing import Tuple, Any, Union

georeader/readers/S2_SAFE_reader.py

@@ -9,13 +9,12 @@
 * Windowed read and read and reproject in the same function (see `load_bands_bbox`)
 * Creation of the image only involves reading one metadata file (`xxx.SAFE/MTD_{self.producttype}.xml`)
 * Compatible with `georeader.read` functions
-* It reads from pyramid if possible
+* It can read from the pyramid if available.
 
 
 https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library
 
 """
-from rasterio import windows
 from shapely.geometry import Polygon
 import xml.etree.ElementTree as ET
 import datetime
@@ -1067,19 +1066,26 @@ def DN_to_radiance(s2obj: S2ImageL1C, dn_data:Optional[GeoTensor]=None) -> GeoTe
      by default this is applied in the S2Image class.
 
      ToA formula from ESA:
-     https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/level-1c/algorithm
+     https://sentiwiki.copernicus.eu/web/s2-processing#S2Processing-TOAReflectanceComputation
 
      Here they say U should be in the numerator
      https://gis.stackexchange.com/questions/285996/convert-sentinel-2-1c-product-from-reflectance-to-radiance
 
      toaBandX = dn_dataBandX / 10000
      radianceBandX = (toaBandX * cos(SZA) * solarIrradianceBandX * U) / pi
 
-     U should be:
-        U = 1. / (1-0.01673*cos(0.0172*(t-4)))^2
-
+     U is:
+        U = 1. / d^2
+    
+    Where d is the Earth-Sun distance correction factor given by:
+        d = 1-0.01673*cos(0.0172*(t-4))
+    
      0.0172 = 360/365.256363 * np.pi/180.
      t = datenum(Y,M,D) - datenum(Y,1,1) + 1;
+     0.01673 is the Earth eccentricity
+
+    In this function we take U from the metadata of the Sentinel-2 image. In the `reflectance` module there are functions
+    to obtain this factor from the date of acquisition of the image.
 
     Args:
         s2obj: s2obj where data has been read

georeader/readers/ee_image.py

@@ -173,6 +173,7 @@ def process_bound(poly):
 def export_cube(query:gpd.GeoDataFrame, geometry:Union[Polygon, MultiPolygon], 
                 transform:Optional[Affine]=None, crs:Optional[str]=None,
                 dtype_dst:Optional[str]=None,
+                bands_gee:Optional[List[str]]=None,
                 crs_polygon:str="EPSG:4326",
                 display_progress:bool=True) -> Optional[GeoTensor]:
     """
@@ -188,6 +189,7 @@ def export_cube(query:gpd.GeoDataFrame, geometry:Union[Polygon, MultiPolygon],
             Defaults to None.
         crs (Optional[str], optional): crs of the geometry. If None it will use the crs of the first image. Defaults to None.
         dtype_dst (Optional[str], optional): dtype of the output GeoTensor. Defaults to None.
+        bands_gee (Optional[List[str]], optional): List of bands to export. If None it will use the bands_gee column in the query. Defaults to None.
         crs_polygon (_type_, optional): crs of the geometry. Defaults to "EPSG:4326".
         display_progress (bool, optional): Display progress bar. Defaults to False.
 
@@ -200,7 +202,9 @@ def export_cube(query:gpd.GeoDataFrame, geometry:Union[Polygon, MultiPolygon],
         return None
 
     # Check required columns
-    required_columns = ["gee_id", "collection_name", "bands_gee"]
+    required_columns = ["gee_id", "collection_name"]
+    if bands_gee is None:
+        required_columns.append("bands_gee")
     if not all([col in query.columns for col in required_columns]):
         raise ValueError(f"Columns {required_columns} are required in the query dataframe")
 
@@ -231,11 +235,17 @@ def export_cube(query:gpd.GeoDataFrame, geometry:Union[Polygon, MultiPolygon],
     def process_query_image(tuple_row):
         _, image = tuple_row
         asset_id = f'{image["collection_name"]}/{image["gee_id"]}'
-        geotensor = export_image_getpixels(asset_id, geometry, proj, image["bands_gee"], crs_polygon=crs_polygon, dtype_dst=dtype_dst)
+        if bands_gee is None:
+            bands_gee_iter = image["bands_gee"]
+        else:
+            bands_gee_iter = bands_gee
+        geotensor = export_image_getpixels(asset_id, geometry, proj, bands_gee_iter, 
+                                           crs_polygon=crs_polygon, dtype_dst=dtype_dst)
         return geotensor
 
     with ThreadPoolExecutor() as executor:
-        geotensor_list = list(tqdm(executor.map(process_query_image, query.iterrows()), total=query.shape[0], disable=not display_progress))
+        geotensor_list = list(tqdm(executor.map(process_query_image, query.iterrows()), 
+                                   total=query.shape[0], disable=not display_progress))
 
     return concatenate(geotensor_list)
 

georeader/readers/ee_query.py

@@ -91,6 +91,27 @@ def figure_out_collection_landsat(tile:str) -> str:
             return "LANDSAT/LC09/C02/T2_TOA"
         else:
             raise ValueError(f"Tile of Landsat-9 {tile} not recognized")
+    elif tile.startswith("LT05"):
+        if tile.endswith("T1"):
+            return "LANDSAT/LT05/C02/T1_TOA"
+        elif tile.endswith("T2"):
+            return "LANDSAT/LT05/C02/T2_TOA"
+        else:
+            raise ValueError(f"Tile of Landsat-5 {tile} not recognized")
+    elif tile.startswith("LT04"):
+        if tile.endswith("T1"):
+            return "LANDSAT/LT04/C02/T1_TOA"
+        elif tile.endswith("T2"):
+            return "LANDSAT/LT04/C02/T2_TOA"
+        else:
+            raise ValueError(f"Tile of Landsat-4 {tile} not recognized")
+    elif tile.startswith("LE07"):
+        if tile.endswith("T1"):
+            return "LANDSAT/LE07/C02/T1_TOA"
+        elif tile.endswith("T2"):
+            return "LANDSAT/LE07/C02/T2_TOA"
+        else:
+            raise ValueError(f"Tile of Landsat-7 {tile} not recognized")
     else:
         raise ValueError(f"Tile {tile} not recognized")
 
@@ -255,6 +276,119 @@ def query(area:Union[MultiPolygon,Polygon],
 
     return geodf
 
+def query_landsat_457(area:Union[MultiPolygon,Polygon],
+                      date_start:datetime, date_end:datetime,
+                      producttype:str="all",
+                      filter_duplicates:bool=True,
+                      return_collection:bool=False,
+                      extra_metadata_keys:Optional[List[str]]=None
+          )-> Union[gpd.GeoDataFrame, Tuple[gpd.GeoDataFrame, ee.ImageCollection]]:
+    """
+    Query Landsat-7, Landsat-5 or Landsat-4 products from the Google Earth Engine.
+
+    Args:
+        area (Union[MultiPolygon,Polygon]): area to query images in EPSG:4326
+        date_start (datetime): datetime in a given timezone. If tz not provided UTC will be assumed.
+        date_end (datetime): datetime in UTC. If tz not provided UTC will be assumed.
+        producttype (str, optional): 'all' -> {"L4", "L5", "L7"}, "L4", "L5" or "L7". Defaults to "all".
+        filter_duplicates (bool, optional): filter duplicate images over the same area. Defaults to True.
+        return_collection (bool, optional): returns also the corresponding image collection. Defaults to False.
+        extra_metadata_keys (Optional[List[str]], optional): extra metadata keys to add to the geodataframe. Defaults to None.
+
+    Returns:
+        Union[gpd.GeoDataFrame, Tuple[gpd.GeoDataFrame, ee.ImageCollection]]: geodataframe with available products in the given area and time range
+    """
+
+    pol = ee.Geometry(mapping(area))
+
+    if date_start.tzinfo is not None:
+        tz = date_start.tzinfo
+        if isinstance(tz, ZoneInfo):
+            tz = tz.key
+
+        date_start = date_start.astimezone(timezone.utc)
+        date_end = date_end.astimezone(timezone.utc)
+    else:
+        tz = timezone.utc
+
+    assert date_end >= date_start, f"Date end: {date_end} prior to date start: {date_start}"
+
+    if extra_metadata_keys is None:
+        extra_metadata_keys = []
+
+    if producttype == "all" or producttype == "L5":
+        image_collection_name = "LANDSAT/LT05/C02/T1_TOA"
+    elif producttype == "L4":
+        image_collection_name = "LANDSAT/LT04/C02/T1_TOA"
+    elif producttype == "L7":
+        image_collection_name = "LANDSAT/LE07/C02/T1_TOA"
+    else:
+        raise NotImplementedError(f"Unknown product type {producttype}")
+
+    keys_query = {"LANDSAT_PRODUCT_ID": "title", 'CLOUD_COVER': "cloudcoverpercentage"}
+    img_col = ee.ImageCollection(image_collection_name).filterDate(date_start.replace(tzinfo=None),
+                                                                   date_end.replace(tzinfo=None)).filterBounds(
+        pol)
+    # Merge T2 collection
+    img_col_t2 = ee.ImageCollection(image_collection_name.replace("T1", "T2")).filterDate(date_start.replace(tzinfo=None),
+                                                                   date_end.replace(tzinfo=None)).filterBounds(
+        pol)
+    img_col = img_col.merge(img_col_t2)
+
+    if producttype == "all":
+        img_col_l4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_TOA").filterDate(date_start.replace(tzinfo=None),
+                                                                   date_end.replace(tzinfo=None)).filterBounds(
+        pol)
+        img_col_l4_t2 = ee.ImageCollection("LANDSAT/LT04/C02/T2_TOA").filterDate(date_start.replace(tzinfo=None),
+                                                                     date_end.replace(tzinfo=None)).filterBounds(
+            pol)
+        img_col_l4 = img_col_l4.merge(img_col_l4_t2)
+        img_col = img_col.merge(img_col_l4)
+
+        # Add L7 T1 and T2
+        img_col_l7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_TOA").filterDate(date_start.replace(tzinfo=None),
+                                                                   date_end.replace(tzinfo=None)).filterBounds(
+        pol)
+        img_col_l7_t2 = ee.ImageCollection("LANDSAT/LE07/C02/T2_TOA").filterDate(date_start.replace(tzinfo=None),
+                                                                        date_end.replace(tzinfo=None)).filterBounds(
+                pol)
+        img_col_l7 = img_col_l7.merge(img_col_l7_t2)
+        img_col = img_col.merge(img_col_l7)
+
+    geodf = img_collection_to_feature_collection(img_col,
+                                                 ["system:time_start"] + list(keys_query.keys()) + extra_metadata_keys,
+                                                as_geopandas=True, band_crs="B2")
+
+    geodf.rename(keys_query, axis=1, inplace=True)
+
+    if geodf.shape[0] == 0:
+        warnings.warn(f"Not images found of collection {producttype} between dates {date_start} and {date_end}")
+        if return_collection:
+            return geodf, img_col
+        return geodf
+
+    img_col = img_col.map(lambda x: _rename_add_properties(x, keys_query))
+    geodf["collection_name"] = geodf.title.apply(lambda x: figure_out_collection_landsat(x))
+
+    geodf = _add_stuff(geodf, area, tz)
+
+    # Fix ids of Landsat to remove initial shit in the names
+    geodf["gee_id"] = geodf["gee_id"].apply(lambda x: "L"+x.split("L")[1])
+
+    if filter_duplicates:
+        geodf = query_utils.filter_products_overlap(area, geodf,
+                                                    groupkey=["solarday", "satellite"]).copy()
+        # filter img_col:
+        img_col = img_col.filter(ee.Filter.inList("title", ee.List(geodf.index.tolist())))
+
+    geodf.sort_values("utcdatetime")
+    img_col = img_col.sort("system:time_start")
+
+    if return_collection:
+        return geodf, img_col
+
+    return geodf
+
 
 def images_by_query_grid(images_available_gee:gpd.GeoDataFrame, grid:gpd.GeoDataFrame) -> gpd.GeoDataFrame:
     aois_images = images_available_gee.reset_index().sjoin(grid, how="inner").reset_index(drop=True)