SAR example¶
Let’s use EOReader with SAR data
Warning:
import os
import glob
# First of all, we need some VHR data, let's use some COSMO-SkyMed data
path = os.path.abspath(glob.glob(os.path.join("/", "*", "DATA", "PRODS", "COSMO", "1011117-766193"))[0])
# Create logger
import logging
logger = logging.getLogger("eoreader")
logger.setLevel(logging.INFO)
# create console handler and set level to debug
ch = logging.StreamHandler()
ch.setLevel(logging.INFO)
# create formatter
formatter = logging.Formatter('%(message)s')
# add formatter to ch
ch.setFormatter(formatter)
# add ch to logger
logger.addHandler(ch)
from eoreader.reader import Reader
# Create the reader
eoreader = Reader()
from eoreader.bands.alias import *
# Open your product
prod = eoreader.open(path)
print(f"Acquisition datetime: {prod.datetime}")
print(f"Condensed name: {prod.condensed_name}")
# Please be aware that EOReader will orthorectify your SAR data with SNAP
# Be sure to have your GPT executable in your path
Acquisition datetime: 2020-10-28 22:46:25
Condensed name: 20201028T224625_CSK_HI_SCS
# Open here some more interesting geographical data: extent
extent = prod.extent()
extent.geometry.to_crs("EPSG:4326").iat[0] # Display
# Open here some more interesting geographical data: footprint
footprint = prod.footprint()
footprint.geometry.to_crs("EPSG:4326").iat[0] # Display
For SAR data, the footprint needs the orthorectified data !
For that, SNAP uses its own DEM, but you can change it when positionning the EOREADER_SNAP_DEM_NAME
environment variable.
Available DEMs are:
ACE2_5Min
ACE30
ASTER 1sec GDEM
Copernicus 30m Global DEM
(buggy for now, do not use it)Copernicus 90m Global DEM
(buggy for now, do not use it)GETASSE30
(by default)SRTM 1Sec HGT
SRTM 3Sec
External DEM
Warning:
If External DEM
is set, you must specify the DEM you want by positioning the EOREADER_DEM_PATH
to a DEM that can be read by SNAP.
# Set the DEM
from eoreader.env_vars import DEM_PATH
os.environ[DEM_PATH] = os.path.join("X:", "BASES_DE_DONNEES", "GLOBAL", "SRTM_30m_v4", "index.vrt")
# Select the bands you want to load
bands = [VV, HH, VV_DSPK, HH_DSPK, HILLSHADE, SLOPE]
# Be sure they exist for COSMO-SkyMed sensor:
ok_bands = [band for band in bands if prod.has_band(band)]
print(to_str(ok_bands)) # This product does not have VV band and HILLSHADE band cannot be computed from SAR band
['HH', 'HH_DSPK', 'SLOPE']
# Load those bands as a dict of xarray.DataArray
band_dict = prod.load(ok_bands)
band_dict[HH]
<xarray.DataArray '20201028T224625_CSK_HI_SCS_HH' (band: 1, y: 16228, x: 16418)> array([[[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32) Coordinates: * band (band) int32 1 * x (x) float64 1.813e+05 1.813e+05 ... 2.306e+05 2.306e+05 * y (y) float64 1.776e+06 1.776e+06 ... 1.727e+06 1.727e+06 spatial_ref int32 0 Attributes: scale_factor: 1.0 add_offset: 0.0 long_name: Intensity_db_HH
xarray.DataArray
'20201028T224625_CSK_HI_SCS_HH'
- band: 1
- y: 16228
- x: 16418
- nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan
array([[[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
- band(band)int321
array([1])
- x(x)float641.813e+05 1.813e+05 ... 2.306e+05
array([181331.469649, 181334.469649, 181337.469649, ..., 230576.469649, 230579.469649, 230582.469649])
- y(y)float641.776e+06 1.776e+06 ... 1.727e+06
array([1776092.033742, 1776089.033742, 1776086.033742, ..., 1727417.033742, 1727414.033742, 1727411.033742])
- spatial_ref()int320
- crs_wkt :
- PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- WGS 84
- horizontal_datum_name :
- World Geodetic System 1984
- projected_crs_name :
- WGS 84 / UTM zone 49N
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- 111.0
- false_easting :
- 500000.0
- false_northing :
- 0.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
- GeoTransform :
- 181329.9696487022 3.0 0.0 1776093.533742319 0.0 -3.0
array(0)
- scale_factor :
- 1.0
- add_offset :
- 0.0
- long_name :
- Intensity_db_HH
# Plot a subsampled version
band_dict[SLOPE][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x1d1006fc898>
# You can also stack those bands
stack = prod.stack(ok_bands)
stack
<xarray.DataArray 'HH_HH_DSPK_SLOPE' (z: 3, y: 16228, x: 16418)> array([[[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ... [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[13.745532, 22.274765, 21.754597, ..., nan, nan, nan], [21.186537, 20.42135 , 19.85021 , ..., nan, nan, nan], [19.413914, 18.580307, 17.949993, ..., nan, nan, nan], ..., [ 0. , 0. , 0. , ..., nan, nan, nan], [ 0. , 0. , 0. , ..., nan, nan, nan], [ 0. , 0. , 0. , ..., nan, nan, nan]]], dtype=float32) Coordinates: spatial_ref int32 0 * x (x) float64 1.813e+05 1.813e+05 ... 2.306e+05 2.306e+05 * y (y) float64 1.776e+06 1.776e+06 ... 1.727e+06 1.727e+06 * z (z) MultiIndex - variable (z) object 'HH' 'HH_DSPK' 'SLOPE' - band (z) int64 1 1 1 Attributes: long_name: ['HH', 'HH_DSPK', 'SLOPE']
xarray.DataArray
'HH_HH_DSPK_SLOPE'
- z: 3
- y: 16228
- x: 16418
- nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan
array([[[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ... [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[13.745532, 22.274765, 21.754597, ..., nan, nan, nan], [21.186537, 20.42135 , 19.85021 , ..., nan, nan, nan], [19.413914, 18.580307, 17.949993, ..., nan, nan, nan], ..., [ 0. , 0. , 0. , ..., nan, nan, nan], [ 0. , 0. , 0. , ..., nan, nan, nan], [ 0. , 0. , 0. , ..., nan, nan, nan]]], dtype=float32)
- spatial_ref()int320
- crs_wkt :
- PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- WGS 84
- horizontal_datum_name :
- World Geodetic System 1984
- projected_crs_name :
- WGS 84 / UTM zone 49N
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- 111.0
- false_easting :
- 500000.0
- false_northing :
- 0.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
- GeoTransform :
- 181329.9696487022 3.0 0.0 1776093.533742319 0.0 -3.0
array(0)
- x(x)float641.813e+05 1.813e+05 ... 2.306e+05
array([181331.469649, 181334.469649, 181337.469649, ..., 230576.469649, 230579.469649, 230582.469649])
- y(y)float641.776e+06 1.776e+06 ... 1.727e+06
array([1776092.033742, 1776089.033742, 1776086.033742, ..., 1727417.033742, 1727414.033742, 1727411.033742])
- z(z)MultiIndex(variable, band)
array([('HH', 1), ('HH_DSPK', 1), ('SLOPE', 1)], dtype=object)
- variable(z)object'HH' 'HH_DSPK' 'SLOPE'
array(['HH', 'HH_DSPK', 'SLOPE'], dtype=object)
- band(z)int641 1 1
array([1, 1, 1], dtype=int64)
- long_name :
- ['HH', 'HH_DSPK', 'SLOPE']
# Error in plotting with a list
if "long_name" in stack.attrs:
stack.attrs.pop("long_name")
# Plot a subsampled version
import matplotlib.pyplot as plt
nrows = len(stack)
fig, axes = plt.subplots(nrows=nrows, figsize=(3*nrows, 6*nrows), subplot_kw={"box_aspect": 1})
for i in range(nrows):
stack[i, ::10, ::10].plot(x="x", y="y", ax=axes[i])