VHR example

Let’s use EOReader with Very High Resolution data.

Warning:

  • We do not provide Pleiades data
  • You will need matplotlib to complete this tutorial
  • import os
    import glob
    
    # First of all, we need some VHR data, let's use Pleiades data
    path = glob.glob(os.path.join("/home", "data", "DATA", "PRODS", "PLEIADES", "5547047101", "IMG_PHR1A_PMS_001"))[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()
    
    # Open your product
    prod = eoreader.open(path, remove_tmp=True)
    print(f"Acquisition datetime: {prod.datetime}")
    print(f"Condensed name: {prod.condensed_name}")
    
    # Please be aware that EOReader will always work in UTM projection, so if you give WGS84 data,
    # EOReader will reproject the stacks and this can be time consuming
    
    Acquisition datetime: 2020-05-11 02:31:58
    Condensed name: 20200511T023158_PLD_ORT_PMS
    
    from eoreader.bands.alias import *
    from eoreader.env_vars import DEM_PATH
    
    # Here, if you want to orthorectify or pansharpen your data manually, you can set your stack here.
    # If you do not provide this stack but you give a non-orthorectified product to EOReader 
    # (ie. SEN or PRJ products for Pleiades), you must provide a DEM to orthorectify correctly the data
    # prod.ortho_stack = ""
    os.environ[DEM_PATH] = os.path.join("/home", "data", "DS2", "BASES_DE_DONNEES", "GLOBAL", "MERIT_Hydrologically_Adjusted_Elevations", "MERIT_DEM.vrt")
    
    # Open here some more interesting geographical data: extent
    extent = prod.extent
    extent.geometry.to_crs("EPSG:4326").iat[0]  # Display
    
    ../_images/VHR_6_01.svg
    # Open here some more interesting geographical data: footprint
    footprint = prod.footprint
    footprint.geometry.to_crs("EPSG:4326").iat[0]  # Display
    
    ../_images/VHR_7_01.svg
    # Select the bands you want to load
    bands = [GREEN, NDVI, TIR_1, CLOUDS, HILLSHADE]
    
    # Be sure they exist for Pleiades sensor:
    ok_bands = [band for band in bands if prod.has_band(band)]
    print(to_str(ok_bands)) # Pleiades doesn't provide TIR and SHADOWS bands
    
    ['GREEN', 'NDVI', 'CLOUDS', 'HILLSHADE']
    
    # Load those bands as a dict of xarray.DataArray
    band_dict = prod.load(ok_bands)
    band_dict[GREEN]
    
    2021-06-09 14:11:49,305 - eoreader - DEBUG - Loading bands ['NIR', 'GREEN', 'RED']
    2021-06-09 14:11:49,305 - eoreader - DEBUG - Read NIR
    2021-06-09 14:11:49,414 - eoreader - INFO - Reprojecting band NIR to UTM with a 0.5 m resolution.
    2021-06-09 14:12:16,191 - eoreader - DEBUG - Manage invalid pixels for band NIR
    2021-06-09 14:12:25,649 - eoreader - DEBUG - Read GREEN
    2021-06-09 14:12:25,760 - eoreader - INFO - Reprojecting band GREEN to UTM with a 0.5 m resolution.
    2021-06-09 14:12:53,321 - eoreader - DEBUG - Manage invalid pixels for band GREEN
    2021-06-09 14:13:02,613 - eoreader - DEBUG - Read RED
    2021-06-09 14:13:02,714 - eoreader - INFO - Reprojecting band RED to UTM with a 0.5 m resolution.
    2021-06-09 14:13:29,244 - eoreader - DEBUG - Manage invalid pixels for band RED
    2021-06-09 14:13:38,462 - eoreader - DEBUG - Loading index ['NDVI']
    2021-06-09 14:13:39,407 - eoreader - DEBUG - Loading DEM bands ['HILLSHADE']
    2021-06-09 14:13:39,407 - eoreader - DEBUG - Warping DEM for IMG_PHR1A_PMS_001
    2021-06-09 14:13:39,409 - eoreader - DEBUG - Using DEM: D:\_EXTRACTEO\DS2\BASES_DE_DONNEES\GLOBAL\MERIT_Hydrologically_Adjusted_Elevations\MERIT_DEM.vrt
    2021-06-09 14:13:43,182 - eoreader - DEBUG - Computing hillshade DEM for IMG_PHR1A_PMS_001
    2021-06-09 14:13:54,877 - eoreader - DEBUG - Loading Cloud bands ['CLOUDS']
    
    <xarray.DataArray 'DIM_PHR1A_PMS_202005110231585_ORT_5547047101' (band: 1, y: 18124, x: 16754)>
    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:
      * x            (x) float64 7.024e+05 7.024e+05 ... 7.108e+05 7.108e+05
      * y            (y) float64 9.689e+06 9.689e+06 9.689e+06 ... 9.68e+06 9.68e+06
      * band         (band) int32 2
        spatial_ref  int32 0
    Attributes: (12/19)
        RADIANCE_BIAS:                         0
        RADIANCE_CALIBRATION_DATE:             2020-04-01T16:00:00.000Z
        RADIANCE_GAIN:                         8.92
        RADIANCE_MEASURE_DESC:                 Raw radiometric count (DN) to TOA ...
        RADIANCE_MEASURE_UNCERTAINTY:          Specification accuracy value
        RADIANCE_MEASURE_UNIT:                 watt/m2/steradians/micrometers
        ...                                    ...
        SPECTRAL_RANGE_MEASURE_DESC:           Spectral range value of raw radiom...
        SPECTRAL_RANGE_MEASURE_UNCERTAINTY:    Specification accuracy value
        SPECTRAL_RANGE_MEASURE_UNIT:           micrometers
        SPECTRAL_RANGE_MIN:                    0.43
        scale_factor:                          1.0
        add_offset:                            0.0
    # The nan corresponds to the nodata you see on the footprint
    %matplotlib inline
    
    # Plot a subsampled version
    band_dict[GREEN][:, ::10, ::10].plot()
    
    <matplotlib.collections.QuadMesh at 0x13a3231ee10>
    
    ../_images/VHR_10_11.png
    # Plot a subsampled version
    band_dict[NDVI][:, ::10, ::10].plot()
    
    <matplotlib.collections.QuadMesh at 0x13a3249e668>
    
    ../_images/VHR_11_11.png
    # Plot a subsampled version
    band_dict[CLOUDS][:, ::10, ::10].plot()
    
    <matplotlib.collections.QuadMesh at 0x13a8e412470>
    
    ../_images/VHR_12_11.png
    # Plot a subsampled version
    band_dict[HILLSHADE][:, ::10, ::10].plot()
    
    <matplotlib.collections.QuadMesh at 0x13a8e4d59b0>
    
    ../_images/VHR_13_11.png
    # You can also stack those bands
    stack = prod.stack(ok_bands)
    stack
    
    2021-06-09 14:14:08,764 - eoreader - DEBUG - Loading bands ['NIR', 'GREEN', 'RED']
    2021-06-09 14:14:08,764 - eoreader - DEBUG - Read NIR
    2021-06-09 14:14:12,928 - eoreader - DEBUG - Read GREEN
    2021-06-09 14:14:16,838 - eoreader - DEBUG - Read RED
    2021-06-09 14:14:20,723 - eoreader - DEBUG - Loading index ['NDVI']
    2021-06-09 14:14:21,655 - eoreader - DEBUG - Loading DEM bands ['HILLSHADE']
    2021-06-09 14:14:21,655 - eoreader - DEBUG - Already existing DEM for IMG_PHR1A_PMS_001. Skipping process.
    2021-06-09 14:14:21,736 - eoreader - DEBUG - Already existing hillshade DEM for IMG_PHR1A_PMS_001. Skipping process.
    2021-06-09 14:14:23,582 - eoreader - DEBUG - Loading Cloud bands ['CLOUDS']
    
    <xarray.DataArray 'NDVI_GREEN_HILLSHADE_CLOUDS' (z: 4, y: 18124, x: 16754)>
    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]],
    
           [[217., 203., 203., ..., 193., 193., 213.],
            [203., 203., 203., ..., 193., 193., 193.],
            [203., 203., 203., ..., 193., 193., 193.],
            ...,
            [ nan,  nan,  nan, ..., 174., 174., 174.],
            [ nan,  nan,  nan, ..., 176., 176., 176.],
            [ nan,  nan,  nan, ..., 177., 177., 196.]],
    
           [[ 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:
        spatial_ref  int32 0
      * x            (x) float64 7.024e+05 7.024e+05 ... 7.108e+05 7.108e+05
      * y            (y) float64 9.689e+06 9.689e+06 9.689e+06 ... 9.68e+06 9.68e+06
      * z            (z) MultiIndex
      - variable     (z) object 'NDVI' 'GREEN' 'HILLSHADE' 'CLOUDS'
      - band         (z) int64 1 1 1 1
    Attributes:
        long_name:  ['NDVI', 'GREEN', 'HILLSHADE', 'CLOUDS']
    # 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=(2*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])
    
    ../_images/VHR_15_01.png