Commit 62c9f538 authored by Julien Michel's avatar Julien Michel
Browse files

Adding Landsat8 dataset and fixing bugs in ecostress dataset and utils

parent cabf8e4e
......@@ -17,12 +17,13 @@ class Ecostress():
"""
ECostress dataset
"""
def __init__(self, lst_file: str, geom_file: str):
def __init__(self, lst_file: str, geom_file: str, cloud_file: str = None):
"""
"""
self.lst_file = lst_file
self.geom_file = geom_file
self.cloud_file = cloud_file
with h5py.File(self.geom_file) as ds:
# Parse acquisition times
......@@ -38,13 +39,13 @@ class Ecostress():
end_time[1:-2])
# Parse bounds
min_lat = ds['StandardMetadata/WestBoundingCoordinate'][()]
max_lat = ds['StandardMetadata/EastBoundingCoordinate'][()]
min_lon = ds['StandardMetadata/SouthBoundingCoordinate'][()]
max_lon = ds['StandardMetadata/NorthBoundingCoordinate'][()]
min_lon = ds['StandardMetadata/WestBoundingCoordinate'][()]
max_lon = ds['StandardMetadata/EastBoundingCoordinate'][()]
min_lat = ds['StandardMetadata/SouthBoundingCoordinate'][()]
max_lat = ds['StandardMetadata/NorthBoundingCoordinate'][()]
self.bounds = rio.coords.BoundingBox(min_lat, min_lon, max_lat,
max_lon)
self.bounds = rio.coords.BoundingBox(min_lon, min_lat, max_lon,
max_lat)
self.crs = '+proj=latlon'
def __repr__(self):
......@@ -154,13 +155,22 @@ class Ecostress():
em[em == 0] = np.nan
vois.append(em)
# Read cloud mask if available
if self.cloud_file:
with h5py.File(self.cloud_file) as cloudDS:
cld = np.array(cloudDS['SDS/CloudMask'][
region[0]:region[2], region[1]:region[3]].astype(dtype))
# CAUTION: we can resample cloud mask with other
# variables as long as we do nearest neighbor
# interpolation
vois.append(cld)
# Stack variables of intereset into a single array
vois = np.stack(vois, axis=-1)
nb_rows = int(np.floor((bounds[2] - bounds[0]) / resolution))
nb_cols = int(np.floor((bounds[3] - bounds[1]) / resolution))
nb_rows = int(np.floor((bounds[3] - bounds[1]) / resolution))
nb_cols = int(np.floor((bounds[2] - bounds[0]) / resolution))
print(nb_rows, nb_cols)
#print(nb_rows, nb_cols)
area_def = pyresample.geometry.AreaDefinition('test', 'test', crs, crs,
nb_cols, nb_rows, bounds)
......@@ -178,16 +188,35 @@ class Ecostress():
lst_end = angles_end + (1 if read_lst else 0)
em_end = lst_end + (5 if read_emissivities else 0)
lst = result[:, :, angles_end] if read_lst else 0
angles = result[:, :, :angles_end] if read_angles else 0
emissivities = result[:, :, lst_end:] if read_emissivities else 0
lst = result[:, :, angles_end] if read_lst else None
angles = result[:, :, :angles_end] if read_angles else None
emissivities = result[:, :, lst_end:] if read_emissivities else None
clouds = result[:, :,
em_end].astype(np.uint8) if self.cloud_file else None
# Unpack cloud mask
masks = None
if self.cloud_file:
valid_mask = np.bitwise_and(clouds, 0b00000001) > 0
cloud_mask = np.logical_or(
np.logical_or(
np.bitwise_and(clouds, 0b00000010) > 0,
np.bitwise_and(clouds, 0b00000100) > 0),
np.bitwise_and(clouds, 0b00001000) > 0)
land_mask = (np.bitwise_and(clouds, 0b00100000) > 0)
sea_mask = np.logical_not(land_mask)
cloud_mask[~valid_mask] = False
land_mask[~valid_mask] = False
sea_mask[~valid_mask] = False
masks = np.stack((cloud_mask, land_mask, sea_mask), axis=-1)
xcoords = np.arange(bounds[0], bounds[0] + nb_cols * resolution,
resolution)
ycoords = np.arange(bounds[3], bounds[3] - nb_rows * resolution,
-resolution)
return lst, emissivities, angles, xcoords, ycoords, crs
return lst, emissivities, angles, masks, xcoords, ycoords, crs
def read_as_xarray(self,
crs: str = None,
......@@ -211,7 +240,7 @@ class Ecostress():
:param dtype: dtype of the output Tensor
"""
lst, emissivities, angles, xcoords, ycoords, crs = self.read_as_numpy(
lst, emissivities, angles, masks, xcoords, ycoords, crs = self.read_as_numpy(
crs, resolution, region, no_data_value, read_lst, read_angles,
read_emissivities, bounds, nprocs, dtype)
......@@ -231,6 +260,11 @@ class Ecostress():
vars['View_Azimuth'] = (['y', 'x'], angles[:, :, 2])
vars['View_Zenith'] = (['y', 'x'], angles[:, :, 3])
if masks is not None:
vars['Cloud_Mask'] = (['y', 'x'], masks[:, :, 0])
vars['Land_Mask'] = (['y', 'x'], masks[:, :, 1])
vars['Sea_Mask'] = (['y', 'x'], masks[:, :, 2])
xarr = xr.Dataset(vars,
coords={
'x': xcoords,
......
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright: (c) 2022 CESBIO / Centre National d'Etudes Spatiales
import os, glob
from typing import List, Union, Tuple
from enum import Enum
import dateutil
import rasterio as rio
import numpy as np
import xarray as xr
from sensorsio import utils
class Landsat:
"""
Class for Landsat L2 product reading
"""
def __init__(self, product_dir: str):
"""
Constructor
:param product_dir: Path to product directory
"""
self.product_dir = os.path.normpath(product_dir)
self.product_name = os.path.basename(self.product_dir)
self.date = dateutil.parser.parse(self.product_name[17:25])
self.year = self.date.year
self.day_of_year = self.date.timetuple().tm_yday
with rio.open(self.build_band_path(Landsat.B1)) as ds:
# Get bounds
self.bounds = ds.bounds
self.transform = ds.transform
# Get crs
self.crs = ds.crs
def __repr__(self):
return f'Landsat {self.date} {self.crs}'
# Enum class for Sentinel2 bands
class Band(Enum):
B1 = 'SR_B1'
B2 = 'SR_B2'
B3 = 'SR_B3'
B4 = 'SR_B4'
B5 = 'SR_B5'
B6 = 'SR_B6'
B7 = 'SR_B7'
B10 = 'ST_B10'
ST_QA = 'ST_QA'
ST_TRAD = 'ST_TRAD'
ST_URAD = 'ST_URAD'
ST_DRAD = 'ST_DRAD'
ST_ATRAN = 'ST_ATRAN'
ST_EMIS = 'ST_EMIS'
ST_EMISD = 'ST_EMISD'
ST_CDIST = 'ST_CDIST'
# Aliases
B1 = Band.B2
B2 = Band.B2
B3 = Band.B3
B4 = Band.B4
B5 = Band.B5
B6 = Band.B6
B7 = Band.B7
B10 = Band.B10
ST_QA = Band.ST_QA
ST_TRAD = Band.ST_TRAD
ST_URAD = Band.ST_URAD
ST_DRAD = Band.ST_DRAD
ST_ATRAN = Band.ST_ATRAN
ST_EMIS = Band.ST_EMIS
ST_EMISD = Band.ST_EMISD
ST_CDIST = Band.ST_CDIST
GROUP_SR = [B1, B2, B3, B4, B5, B6, B7]
GROUP_ST = [B10]
class Mask(Enum):
CLOUDS = 'QA_PIXEL'
AEROSOLS = 'QA_AEROSOL'
SATURATIONS = 'QA_RADSAT'
# Aliases
CLOUDS = Mask.CLOUDS
AEROSOLS = Mask.AEROSOLS
SATURATIONS = Mask.SATURATIONS
ALL_MASKS = [CLOUDS, AEROSOLS, SATURATIONS]
# From https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/LSDS-1619_Landsat8-C2-L2-ScienceProductGuide-v2.pdf
FACTORS = {
B1: 0.0000275,
B2: 0.0000275,
B3: 0.0000275,
B4: 0.0000275,
B5: 0.0000275,
B6: 0.0000275,
B7: 0.0000275,
B10: 0.00341802,
ST_QA: 0.01,
ST_TRAD: 0.001,
ST_URAD: 0.001,
ST_DRAD: 0.001,
ST_ATRAN: 0.0001,
ST_EMIS: 0.0001,
ST_EMISD: 0.0001,
ST_CDIST: 0.01
}
SHIFTS = {
B1: -0.2,
B2: -0.2,
B3: -0.2,
B4: -0.2,
B5: -0.2,
B6: -0.2,
B7: -0.2,
B10: 149,
ST_QA: 0,
ST_TRAD: 0,
ST_URAD: 0,
ST_DRAD: 0,
ST_ATRAN: 0,
ST_EMIS: 0,
ST_EMISD: 0,
ST_CDIST: 0
}
NO_DATA_FLAGS = {
B1: 0,
B2: 0,
B3: 0,
B4: 0,
B5: 0,
B6: 0,
B7: 0,
B10: 0,
ST_QA: -9999,
ST_TRAD: -9999,
ST_URAD: -9999,
ST_DRAD: -9999,
ST_ATRAN: -9999,
ST_EMIS: -9999,
ST_EMISD: -9999,
ST_CDIST: -9999
}
def build_band_path(self, band: Union[Band, Mask]) -> str:
"""
Build path to a band for product
:param band: The band to build path for as a Sentinel2.Band enum value
:return: The path to the band file
"""
p = glob.glob(f"{self.product_dir}/*{band.value}.TIF")
# Raise
if len(p) == 0:
raise FileNotFoundError(
f"Could not find band {band.value} in product directory {self.product_dir}"
)
return p[0]
def read_as_numpy(
self,
bands: List[Band],
masks: List[Mask] = ALL_MASKS,
crs: str = None,
resolution: float = 30,
region: Union[Tuple[int, int, int, int],
rio.coords.BoundingBox] = None,
no_data_value: float = np.nan,
bounds: rio.coords.BoundingBox = None,
algorithm=rio.enums.Resampling.cubic,
dtype: np.dtype = np.float32
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, str]:
"""
Read bands from Sentinel2 products as a numpy ndarray. Depending on the parameters, an internal WarpedVRT
dataset might be used.
:param bands: The list of bands to read
:param crs: Projection in which to read the image (will use WarpedVRT)
:param resolution: Resolution of data. If different from the resolution of selected bands, will use WarpedVRT
:param region: The region to read as a BoundingBox object or a list of pixel coords (xmin, ymin, xmax, ymax)
:param no_data_value: How no-data will appear in output ndarray
:param bounds: New bounds for datasets. If different from image bands, will use a WarpedVRT
:param algorithm: The resampling algorithm to be used if WarpedVRT
:param dtype: dtype of the output Tensor
:return: The image pixels as a np.ndarray of shape [bands, width, height],
The x coords as a np.ndarray of shape [width],
the y coords as a np.ndarray of shape [height],
the crs as a string
"""
np_arr = None
np_arr_msk = None
xcoords = None
ycoords = None
crs = None
# Readn bands
if len(bands):
img_files = [self.build_band_path(b) for b in bands]
np_arr, xcoords, ycoords, crs = utils.read_as_numpy(
img_files,
crs=crs,
resolution=resolution,
region=region,
output_no_data_value=no_data_value,
bounds=bounds,
algorithm=algorithm,
separate=True,
dtype=dtype)
factors = np.array([self.FACTORS[b] for b in bands])
shifts = np.array([self.SHIFTS[b] for b in bands])
# Skip first dimension
np_arr = np_arr[0, ...]
np_arr_rescaled = (factors * np_arr) + shifts
for i, b in enumerate(bands):
np_arr_rescaled[i, ...][np_arr[i, ...] ==
self.NO_DATA_FLAGS[b]] = no_data_value
np_arr = np_arr_rescaled
if len(masks):
img_files = [self.build_band_path(m) for m in masks]
np_arr_msk, xcoords, ycoords, crs = utils.read_as_numpy(
img_files,
crs=crs,
resolution=resolution,
region=region,
output_no_data_value=no_data_value,
bounds=bounds,
algorithm=rio.enums.Resampling.nearest,
separate=True,
dtype=np.uint16,
scale=None)
# Drop first dimension
np_arr_msk = np_arr_msk[0, ...]
return np_arr, np_arr_msk, xcoords, ycoords, crs
def read_as_xarray(self,
bands: List[Band],
masks: List[Mask] = ALL_MASKS,
crs: str = None,
resolution: float = 30,
region: Union[Tuple[int, int, int, int],
rio.coords.BoundingBox] = None,
no_data_value: float = np.nan,
bounds: rio.coords.BoundingBox = None,
algorithm=rio.enums.Resampling.cubic,
dtype: np.dtype = np.float32) -> xr.Dataset:
"""
Read bands from Sentinel2 products as a numpy ndarray. Depending on the parameters, an internal WarpedVRT
dataset might be used.
:param bands: The list of bands to read
:param crs: Projection in which to read the image (will use WarpedVRT)
:param resolution: Resolution of data. If different from the resolution of selected bands, will use WarpedVRT
:param region: The region to read as a BoundingBox object or a list of pixel coords (xmin, ymin, xmax, ymax)
:param no_data_value: How no-data will appear in output ndarray
:param bounds: New bounds for datasets. If different from image bands, will use a WarpedVRT
:param algorithm: The resampling algorithm to be used if WarpedVRT
:param dtype: dtype of the output Tensor
:return:
"""
np_arr, np_arr_msk, xcoords, ycoords, crs = self.read_as_numpy(
bands, masks, crs, resolution, region, no_data_value, bounds,
algorithm, dtype)
vars = {}
for i in range(len(bands)):
vars[bands[i].value] = (["t", "y", "x"], np_arr[None, i, ...])
for i in range(len(masks)):
vars[masks[i].value] = (["t", "y", "x"], np_arr_msk[None, i,
...])
xarr = xr.Dataset(vars,
coords={
't': [self.date],
'x': xcoords,
'y': ycoords
},
attrs={'crs': crs})
return xarr
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright: (c) 2021 CESBIO / Centre National d'Etudes Spatiales
"""
This module contains utilities function
"""
......@@ -17,16 +16,14 @@ from rasterio.warp import transform_bounds
import rasterio as rio
from affine import Affine
def rgb_render(data: np.ndarray,
clip: int = 2,
bands: List[int] = [2,
1,
0],
norm: bool = True,
dmin: np.ndarray = None,
dmax: np.ndarray = None) -> Tuple[np.ndarray,
np.ndarray,
np.ndarray]:
def rgb_render(
data: np.ndarray,
clip: int = 2,
bands: List[int] = [2, 1, 0],
norm: bool = True,
dmin: np.ndarray = None,
dmax: np.ndarray = None) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Prepare data for visualization with matplot lib
......@@ -37,8 +34,8 @@ def rgb_render(data: np.ndarray,
:returns: a tuple of data ready for matplotlib, dmin, dmax
"""
assert(len(bands) == 1 or len(bands) == 3)
assert(clip >= 0 and clip <= 100)
assert (len(bands) == 1 or len(bands) == 3)
assert (clip >= 0 and clip <= 100)
# Extract bands from data
data_ready = np.take(data, bands, axis=0)
......@@ -63,11 +60,10 @@ def rgb_render(data: np.ndarray,
return data_ready, dmin, dmax
def generate_psf_kernel(
res: float,
mtf_res: float,
mtf_fc: float,
half_kernel_width: int = None) -> np.ndarray:
def generate_psf_kernel(res: float,
mtf_res: float,
mtf_fc: float,
half_kernel_width: int = None) -> np.ndarray:
"""
Generate a gaussian PSF kernel sampled at given resolution
......@@ -99,16 +95,15 @@ def generate_psf_kernel(
return kernel
def create_warped_vrt(
filename: str,
resolution: float,
dst_bounds: BoundingBox = None,
dst_crs: str = None,
src_nodata: float = None,
nodata: float = None,
shifts: Tuple[float] = None,
resampling: Resampling = Resampling.cubic,
dtype=None) -> WarpedVRT:
def create_warped_vrt(filename: str,
resolution: float,
dst_bounds: BoundingBox = None,
dst_crs: str = None,
src_nodata: float = None,
nodata: float = None,
shifts: Tuple[float] = None,
resampling: Resampling = Resampling.cubic,
dtype=None) -> WarpedVRT:
"""
Create a warped vrt from filename, to change srs and resolution
......@@ -136,7 +131,7 @@ def create_warped_vrt(
target_bounds = transform_bounds(src.crs, dst_crs, *src.bounds)
else:
target_bounds = src.bounds
src_transform = src.transform
if shifts is not None:
src_res = src_transform[0]
......@@ -147,8 +142,7 @@ def create_warped_vrt(
left, bottom, right, top = target_bounds
dst_width = (right - left) / resolution
dst_height = (top - bottom) / resolution
dst_transform = Affine(resolution, 0.0, left,
0.0, -resolution, top)
dst_transform = Affine(resolution, 0.0, left, 0.0, -resolution, top)
vrt_options = {
'resampling': resampling,
......@@ -157,7 +151,6 @@ def create_warped_vrt(
'width': dst_width,
'crs': target_crs,
'src_transform': src_transform
}
if src_nodata is not None:
vrt_options['src_nodata'] = src_nodata
......@@ -173,6 +166,7 @@ def create_warped_vrt(
return vrt
def bb_intersect(bb: List[BoundingBox]) -> BoundingBox:
"""
Compute the intersection of a list of bounding boxes
......@@ -192,6 +186,7 @@ def bb_intersect(bb: List[BoundingBox]) -> BoundingBox:
return BoundingBox(left=xmin, bottom=ymin, right=xmax, top=ymax)
def bb_snap(bb: BoundingBox, align: float = 20) -> BoundingBox:
"""
Snap a bounding box to multiple of align parameter
......@@ -207,7 +202,11 @@ def bb_snap(bb: BoundingBox, align: float = 20) -> BoundingBox:
top = align * np.ceil(bb[3] / align)
return BoundingBox(left=left, bottom=bottom, right=right, top=top)
def bb_common(bounds: List[BoundingBox], src_crs:List[str], snap: float = 20, target_crs: str = None):
def bb_common(bounds: List[BoundingBox],
src_crs: List[str],
snap: float = 20,
target_crs: str = None):
"""
Compute the common bounding box between a set of images.
All bounding boxes are converted to crs before intersection.
......@@ -221,32 +220,34 @@ def bb_common(bounds: List[BoundingBox], src_crs:List[str], snap: float = 20, ta
returns: A tuple of box, crs
"""
assert(len(bounds)==len(src_crs))
boxes=[]
assert (len(bounds) == len(src_crs))
boxes = []
for box, crs in zip(bounds, src_crs):
if target_crs is None:
target_crs=crs
target_crs = crs
crs_box = rio.warp.transform_bounds(crs, target_crs, *box)
boxes.append(crs_box)
# Intersect all boxes
box = bb_intersect(boxes)
# Snap to grid
box = bb_snap(box, align=snap)
return box, crs
def read_as_numpy(img_files:List[str],
crs: str=None,
resolution:float = 10,
offsets:Tuple[float,float]=None,
region:Union[Tuple[int,int,int,int],rio.coords.BoundingBox]=None,
input_no_data_value:float=None,
output_no_data_value:float=np.nan,
bounds:rio.coords.BoundingBox=None,
return box, target_crs
def read_as_numpy(img_files: List[str],
crs: str = None,
resolution: float = 10,
offsets: Tuple[float, float] = None,
region: Union[Tuple[int, int, int, int],
rio.coords.BoundingBox] = None,
input_no_data_value: float = None,
output_no_data_value: float = np.nan,
bounds: rio.coords.BoundingBox = None,
algorithm=rio.enums.Resampling.cubic,
separate:bool=False,
separate: bool = False,
dtype=np.float32,
scale:float=None) -> np.ndarray:
scale: float = None) -> np.ndarray:
"""
:param vrts: A list of WarpedVRT objects to stack
:param region: The region to read as a BoundingBox object or a list of pixel coords (xmin, ymin, xmax, ymax)
......@@ -255,62 +256,70 @@ def read_as_numpy(img_files:List[str],
TODO
"""
"""
#print(f'{bounds=}')
# Check if we need resampling or not
need_warped_vrt = (offsets is not None)
# If we change image bounds
for f in img_files:
with rio.open(f) as ds:
if bounds is not None and ds.bounds != bounds:
need_warped_vrt=True
need_warped_vrt = True
# If we change projection