Orthorectification and radiometric calibration of raw Pleiades data yields unreliable results - What am I doing wrong?

[eikeschuett]

I got access to "raw" Pleiades data (Geometric Processing Level: SENSOR; Radiometric Processing Level: BASIC), which I'm trying to orthorectify, calibrate radiometrically and export as GeoTIFF. I'm new to Pleiades data and was lucky finding EOReader to ease this process. If I understood correctly, opening the product, specifying a DEM and creating a stack should do all this. However, the final GeoTIFF appears to be offset by ~12 m towards ENE, if compared with e.g. Google or Bing Maps.

Additionally, the reflectances of e.g. the Blue band are in a range between 411 and 4267, which appears quite strange to me. What is the unit of this output? Or may this be related to an error in the radiometric calibration? Strangely, if I include the {TO_REFLECTANCE: False} keyword in the prod.stack method, the GeoTIFF contains only white noise.

It may be that both things are results of improper usage of your package. I built my code mainly from snippets from different tutotials:

import logging
import tempenv

from sertit import logs

from eoreader.reader import Reader
from eoreader.bands import GREEN, BLUE, RED, NIR, DEM
from eoreader.env_vars import DEM_PATH
from eoreader.keywords import TO_REFLECTANCE
from dask.distributed import Client, LocalCluster

def run():
    # Create logger
    logger = logging.getLogger("eoreader")
    logs.init_logger(logger)
    
    pleiades_fname = "IMG_PHR1A_PMS_001"
    
    reader = Reader()
    
    prod = reader.open(pleiades_fname)
    
    try:
        prod.load(RED)[RED]

    except ValueError as msg:
       logger.error(msg)
       
    dem = "SH_UTM_WGS84.tif"
    
    with LocalCluster(
            n_workers=4, threads_per_worker=4, processes=True
                      ) as cluster, Client(cluster):
        
        with tempenv.TemporaryEnvironment({DEM_PATH: dem}):
            prod.load(DEM)[DEM]
    
            stack = prod.stack([BLUE], 
                               stack_path="test_stack.tif",
                               save_as_int = True)
            
    return stack
   
if __name__ == "__main__":
    stack = run()

Environment:

• Python version: 3.10.6
• eoreader version: 0.16.0

Further Ideas
I plan to use Pleiades data for aquatic purposes. Therefore, I'm more interested in at-sensor radiances than reflectances to perform atmospheric correction and then later calculate Remote Sensing Reflactances (R_rs) myself. I haven't found a possibility to skip the calculation of reflectances (except editing the source code). Including a keyword for that purpose would be a great feature, I think!

[remi-braun]

Hello,

Thanks for this extra-detailed feedback!
You might have found some bugs 😅

First of all, what kind of DEM is SH_UTM_WGS84 ? What is its resolution ? I never encountered such displacement.
Could you retry using another one ? (maybe COPDEM or SRTM ?)
Side Note: In the provided snippet, if you don't need the DEM band, you don't need to load it 😉

Concerning the reflectance & white noise:

• The shown integar values are those stored on disk or loaded in Python ? If it comes from the disk values, it is normal as the values are saved as uint16. Could you try without save_as_int argument ? I may have a bug with TO_REFLECTANCE = False and save_as_int = True though.
• If this doesn't resolve your ptoblem, it may come from Dask. Could you try without any cluster ? Dask handling is very very preliminar, we (still) don't really use it, and therefore might be a bit buggy...
• If you are still in a deadlock, could you give me what is printed in the logs ? And maybe the DIMAP.xml file to investigate the metadata ? It could come from a Pleiades case I never have seen before.

To answer your last idea, the keyword to keep the bands "as is" is already implemented and is TO_REFLECTANCE = False 😉

[eikeschuett]

Thanks for your answer, Rémi!

• Concerning the reflectance issues: Without the save_as_int argument, everything works fine. Both for DNs and for reflectances. That is great!

• Regarding the spatial offset: I used the EU-DEM with a spatial resolution of 25 m. I tried it again with SRTM and another scene, but there is still a considerable offset. I should point out again that the Geometric Processing Level of the data is "Sensor" (i.e. still in sensor geomety), but I thought that eoreader takes care of the projection and orthorectification, right?
  

• Regarding the output as radiance: I don't want to keep my band "as is", since my data provides raw DN. I want to do the radiometric calibration from DNs to radiances and output the radiances (no further calculation of the reflectances). I had a brief look into dimap_product.py and I guess we would need a new method inside the DimapProductclass, which basically does exactly the same as _to_reflectance(), except line 667. I guess there would be a smarter way to implement this (if you think this option would be useful. For me it would be great!)...

I think eoreader is a great project, which will ease the workflows of my new project a lot. Thanks for your effort and time!

[remi-braun]

Thanks for trying the workarounds!

➡️ The issue when saving stacks as uint16 with bands loaded "as is" is referenced here: #52

➡️ I think I never really paid attention to SEN products (we don't use them during CEMS activations).
What I'm doing is I am reprojecting the stacks with the given RPCs on the DEM.
I also have some shifts on my SEN products, but I was assuming it came from the DEM I'm using (orthorectifying a 0.5m image with a 90m DEM for testing purposes).
But since you have also this issue, I have referenced it here: #53
I have a workaroud though. If you can orthorectify the stack elsewhere, just set the ortho stack path as described here:

Here, if you want to orthorectify or pansharpen your data manually, you can set your stack here.

prod.ortho_stack = "/path/to/ortho_stack.tif"

➡️ Sorry I read too fast for the radiance output. You are right, for now you cannot load bands as radiance in EOReader.
This won't be too difficult to implement, just a bit long as I need to do it for every handled product. And since we don't need it in CEMS activations, this won't be my top priority 😅. The issue is referenced here: #54
The workaround for you could be loading the bands "as is" and then using the following private method (this isn't tested but should work):

# Convert DN into radiance
nir_array_rad = prod._dn_to_toa_rad(nir_array_dn, NIR)

⚠️ This methods doesn't exist for every optical products!

Feel free to answer directly in the corresponding issues 😉