OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.84k stars 2.53k forks source link

.nat format of MeTOP AVHRR L1B not readable #5627

Open Scartography opened 2 years ago

Scartography commented 2 years ago

In the GDAL driver list stands that the .nat format should be readable for NOAA AVHRR is readable, which might be true, however for the similar data source MeTOP AVHRR GDAL (GDAL 3.4.0, released 2021/11/04) cannot read these files:

$ gdalinfo AVHR_xxx_1B_M01_20210710081003Z_20210710095203Z_N_O_20210710085555Z.nat ERROR 4: AVHR_xxx_1B_M01_20210710081003Z_20210710095203Z_N_O_20210710085555Z.nat' not recognized as a supported file format.

The data can be found at: https://data.eumetsat.int/data/map/EO:EUM:DAT:METOP:AVHRRL1

I would have assumed that it would be readable as per: https://gdal.org/drivers/raster/msgn.html

As I have just started dealing with these type of data, I might be wrong on this. Here is a script by which I now search and handle the data as first dev test. This for now is unrelated to GDAL, but might help anyone here to get the data without using the browser.

import shutil
from glob import glob

from satpy import Scene
from satpy import available_readers
from satpy import DataQuery
from datetime import datetime

import eumdac

import os
os.environ['PYTROLL_CHUNK_SIZE'] = "1024"
os.environ['DASK_NUM_WORKERS'] = "2"
os.environ['OMP_NUM_THREADS'] = "1"

consumer_key = "**"*
consumer_secret = "***"

credentials = (consumer_key, consumer_secret)

token = eumdac.AccessToken(credentials)

print(f"This token '{token}' expires {token.expiration}")

datastore = eumdac.DataStore(token)
# print(datastore.collections)

selected_collection = datastore.get_collection('EO:EUM:DAT:METOP:AVHRRL1')

# Add vertices for polygon, wrapping back to the start point.
geometry = [
    [-1.0, -1.0],
    [4.0, -4.0],
    [8.0, -2.0],
    [9.0, 2.0],
    [6.0, 4.0],
    [1.0, 5.0],
    [-1.0, -1.0]
    ]

# Set sensing start and end time
start = datetime(2021, 7, 10, 8, 0)
end = datetime(2021, 7, 10, 9, 0)

# Retrieve datasets that match our filter
products = selected_collection.search(
    geo='POLYGON(({}))'.format(
        ','.join(["{} {}".format(*coord) for coord in geometry])
        ),
    dtstart=start,
    dtend=end
    )

print(f'Found Datasets: {len(products)} datasets for the given time range')

# for product in products:
#     with product.open() as fsrc, \
#             open(fsrc.name, mode='wb') as fdst:
#         shutil.copyfileobj(fsrc, fdst)

# print(available_readers())

satpy_products = []
for product in products:
    satpy_products.append(
        "_path_to_extracted_products"
        f"{str(product)}/{str(product)}.nat"
        )

scn = Scene(reader="avhrr_l1b_eps", filenames=satpy_products)
print(scn.all_dataset_ids())
scn.load([0.58, 0.63, 0.68], calibration=['reflectance'])
new_sc = scn.resample(resampler='native')
new_sc.save_datasets(writer='geotiff')

print(scn.all_dataset_ids())
# print(filenames)
# my_channel_id = DataQuery(name='IR_016', calibration='reflec_1')

# print(scn.load(calibration=["radiance"]))
rouault commented 2 years ago

Looking at one sample MeTOP AVHRR L1B file and the MSGN driver, they appear to be in a different format whose specification in https://www.eumetsat.int/media/38675