Extent not detected from lon & lat vars in netcdf files from Acolite AC products

[arthurdegrandpre]

Hi,

I have been using terra to work with satellite imagery, and have come across issues when reading some netcdf products.

Specifically, when opening outputs from Acolite atmospheric correction where lon and lat are stored as vars, terra does not recognize the extent of the data. In some cases, an additionnal var is included with the Acolite products containing projection information, preventing the issue.

When the first two vars are lon and lat, without any projection var, this results in the following warning message :
"vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named x BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"

And the following SpatRaster properties :
extent : 0.2840993, 17542.72, 0.5, 17542.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs

Neither lon or lat vars are included in the bands of the SpatRaster object.

I am able to bypass this issue by using the ncdf4 library the following way:

image = ncdf4::nc_open(image_path)

lon = ncdf4::ncvar_get(image, image$var[[1]])
lat = ncdf4::ncvar_get(image, image$var[[2]])
b1 = ncdf4::ncvar_get(image, image$var[[3]])

r = terra::rast(b1, ext=c(range(lon), range(lat)), crs="+proj=lonlat")

Is there any way to force terra to read the product into a valid georeferencing? (In this case that would be WGS84)

Many thanks,
Arthur

[mdsumner]

Can you point to specific products, and specific urls we can download or access via streaming (requiring auth is ok). This is very common, and GDAL has a lot of heuristics but not for every actual let alone imaginable case of how these arrays can be structured (sometimes the coordinate(geolocation)-arrays are real, and you must go through the warper to resolve them to a regular grid, and sometimes the geolocation arrays aren't registered to GDAL, and sometimes they are entirely redundant. It's not possible to determine which is the case automatically).

[rhijmans]

I would love to help, but I would need an example file.

What you are doing here:

r = terra::rast(b1, ext=c(range(lon), range(lat)), crs="+proj=lonlat")

Is probably not correct because the extent is not defined by the coordinates of centers of the grid cells.

[arthurdegrandpre]

Hello and thanks to both of you,
First of all sorry for the delayed answer. There is a lot to unpack here and I lacked time for a deep dive.

I now realize this issue might actually be better directed towards GDAL rather than terra, or maybe even @acolite. Sorry if this is misguided. I did suspect an issue with GDAL being unable to parse this specific metadata structure, although I have been able to read the products directly through SNAP which I believe also uses GDAL.

Specifically, in my situation, the original data product are Pleiades Neo products generated by AIRBUS DS GEO, accessed through the Skywatch platform (a third party distributor). I do not have permission to share the products so that makes things harder.
I will verify if I could setup a shared folder containing related metadata files, maybe an example .nc file output from Acolite where reflectance values are wiped. I could also send Acolite log and settings files.

As that might not be possible, here is some information regarding the products and their processing:
The original data source (one of many) is type ORTHO, 6 bands reflectance product from dataset ID ACQ_PNEO3_04734513309190, product name PNEO3_202408291550514_MS-FS_ORT in DIMAP format. Data structure for one product is 2 folders for MS and PAN data, with an INDEX and VOL_PNEO html file at the root, and another set of INDEX, DIM, LUT and RPC HTML files, and series of .TIF and .TFW files for both RBG and NED bands per tile. There are also .GML files in a MASKS folder, and more XML metadata in a LINEAGE folder.

This image was processed using ACOLITE (see https://github.com/acolite/acolite) into L2W water leaving reflectance into .nc files, where lon and lat are stored as vars, but does not seem to contain CRS information. It appears that ACOLITE itself already has difficulty parsing the geolocation information. While this processed product is readable in SNAP and seemingly well localized, geolocation data does not parse properly into terra. I was able to manually unpack the .nc file into a spatRaster using ncdf4 and terra by manually copying the metadata available in the INDEX metadata file from the original Airbus/Skywatch product, but I wasn't able to programmatically parse the INDEX files for the information.

@rhijmans you are right that the previously mentioned method didn't produce accurate results. I have adjusted my scripts accordingly. Here is an update snippet, if it might be of use to anyone, where values were manually copied from original metadata files.

`
nc = ncdf4::nc_open(image_path)
lon = ncdf4::ncvar_get(nc, nc$var[[which(names(nc$var) == "lon")]])

image_nc = terra::rast(nrows=dim(lon)[2],
ncols=dim(lon)[1],
xmin = copy_from_INDEX,
xmax = copy_from_INDEX+res_from_INDEXdim(lon)[1],
ymax = copy_from_INDEX,
ymin = copy_from_INDEX-res_from_INDEXdim(lon)[2],
resolution = res_from_INDEX,
crs = "EPSG_code_from_INDEX")
`

Thanks again

[acolite]

Hi Arthur

The Pléiades support in ACOLITE is experimental (and especially for Neo) as I don't have access to very many images, and there are quite a few configuration options for the DIMAP bundles. My feeling is that an ORTHO dataset should have a projection associated, so I am not exactly sure what is happening with your files. What does gdalinfo output for the .TIF files in your DIMAP bundle?

ACOLITE tries to read the projection information from the input image (using gdal). If this is successful, the projection information is (should be...) stored in the NetCDF. If there is no projection information in the image, the geolocation is estimated by simple linear interpolation from the image corners and centre provided in the metadata. In this case there is no projection information in the NetCDF, but lat and lon datasets are written anyway.

In the ACOLITE processing there will be a message in your log file that reads "Could not determine projection from ..." and the processing will continue with "Pléiades old method". If the projection data could be read then the message will be "Pléiades new method" instead. (Those are likely the files you have no issue with? However I only just now committed the projection info being stored in the NetCDF for the "new method". Do a git pull for best results.)

If there are RPC provided with the image, then you can reproject the input and have that projection information stored in the NetCDF. (Note that image dimensions will change from the input file.) For the files with no associated projection, try the settings:

reproject_inputfile=True 
reproject_inputfile_dem=True

For sea level images, the reproject_inputfile_dem is not really needed and can be omitted. This option retrieves Copernicus DEM data that covers the image and passes this as RPC_DEM in the transformerOptions to the gdal.Warp function.

I hope this helps!

Quinten

[arthurdegrandpre]

Hi Quentin, thanks for the quick feedback.

Running gdalinfo does appear to retrieve the whole of the Coordinate system properly.

GDAL INFO PNEO.txt

In both the newer Acolite version and the one I previously used, the Pléiades new method is used.
The pulled Acolite version with the reproject_inputfile options set to True did result in the correct identification of the CRS when reading from terra and the loss of the warning message:

"vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named y BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"

However, the extent origin is still (0,0) and the resolution is set to 1,1 instead of the expected 1.2,1.2 in the spatRaster.
Examining the .nc file shows a new var called transverse_mercator but neither by examining the vars or attributes was I able to locate the missing information.

Based on these observations I think it might be appropriate to transfer this issue towards the Acolite Issues page.

[acolite]

Did you try the updated version without reproject_inputfile?

Quinten

[arthurdegrandpre]

Changing reproject_inputfile from True to False did not change the outcome.
Also, reading in terra generates the following error message:

terra::rast(image_path)

Warning message:
[rast] cells are not equally spaced; extent is not defined

[acolite]

Do the tif outputs work properly when you run with l2w_export_geotiff=True?

Quinten

[arthurdegrandpre]

The tif outputs are generated and the projection information is stored inside, but the extent is missing so the origin remains 0,0 and the spatial resolution is also set to 1 (both read using terra and QGIS). Acolite log doesn't show any particular errors or warnings.

Thanks,
Arthur

[acolite]

Is this without reproject_inputfile?

In that case the resolution info should be copied from the original tif.

Quinten

[arthurdegrandpre]

Both outputs from reproject_inputfile True and False result in the same extent and resolution issue.

Arthur

[acolite]

What does gdalinfo give for both outputs?

Could you share the original DIMAP bundle with blanked tiffs?

Quinten

[arthurdegrandpre]

Here are the gdalinfo files, I'll work on blanking the tiffs.
Arthur

gdalinfo_reproject_false.txt
gdalinfo_reproject_true.txt

[arthurdegrandpre]

Here is the blanked bundle, where all values were set to 1. I left only a single tile and did not include the PAN data. I made sure Acolite was still able to process the bundle and replicate the issue.

blanked_000248901_1_2_STD_A.zip

Thanks again,
Arthur

[acolite]

Hi Arthur

It seems that QGIS (GDAL) can position the ACOLITE NetCDF and GeoTIFF outputs correctly for this blanked image, see the screenshot below with a NetCDF layer and a GeoTIFF layer overlapping, and WorldLakes in the background. Extent, origin, and pixel size are all populated in the Layer properties window.

Here is the GeoTIFF output if that helps anyone:

PNEO4_2024_09_11_15_50_26_L1R_rhot_826.tif.zip

Quinten