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

[arthurdegrandpre]

Thank you for running the test on your end Quinten. Sadly, the issue persists on my end when processing the same blank product I provided... This issue does not occur with products from other sources, such as planet.

May I have a copy of your run's l1r_settings? Are there any parameters that could cause issues with PNEO only?

Arthur

[acolite]

Did the tif I generated work properly for you?

I set no special parameters. All I had was:

inputfile=blanked_000248901_1_2_STD_A/
output=Output/
l1r_export_geotiff=True
l2r_export_geotiff_rgb=True

Quinten

[arthurdegrandpre]

It does work! It made me realize I never checked the L1R and L2R products, as I they are deleted in favor of L2W. I changed my settings to keep them and they also work fine. Only the L2W products seem to have the extent issue...

Arthur

[mdsumner]

I'm still confident that this is a netcdf-variant issue. NetCDF is not modelled around raster-logic with an efficient and high-fidelity extent specification. Can you point to an actual L2W file? I'll provide a workaround for terra, and report/fix to GDAL as well if anything can be improved. (Many things have been improved in GDAL recently, so it might already be fixed and just need an interim workaround for now). We can't fix the broken world of NetCDF, though.

[acolite]

Here is a sample L2W NetCDF file and L2W GeoTIFF exported from that NetCDF:

PNEO4_2024_09_11_15_50_26_L2W.nc.zip PNEO4_2024_09_11_15_50_26_L2W_rhos_826.tif.zip

Arthur, does this TIF also works for you?

[mdsumner]

I got this to work in a way that suprised me a bit, thanks for the example! (I can't speak on @rhijmans behalf about terra but I can navigate the GDAL relationship in a way that is also of active interest and ongoing related-work for me).

rast("/vsizip//vsicurl/https://github.com/user-attachments/files/20721164/PNEO4_2024_09_11_15_50_26_L2W.nc.zip", "rhos_826", drivers = "NetCDF")
class       : SpatRaster
dimensions  : 6742, 7680, 1  (nrow, ncol, nlyr)
resolution  : 1.2, 1.2  (x, y)
extent      : 657966, 667182, 5108896, 5116986  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
source      : PNEO4_2024_09_11_15_50_26_L2W.nc.zip:rhos_826
varname     : rhos_826 (Surface reflectance)
name        : rhos_826
unit        :        1

Using vsicurl to read from a URL or vsizip from a zip won't work on windows (or macos I think). There we can use this by download the file, as for me on windows (where GDAL is 3.10.2):

rast("PNEO4_2024_09_11_15_50_26_L2W.nc", "rhos_826", drivers = "NetCDF")
class       : SpatRaster 
dimensions  : 6742, 7680, 1  (nrow, ncol, nlyr)
resolution  : 1.2, 1.2  (x, y)
extent      : 657966, 667182, 5108896, 5116986  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
source      : PNEO4_2024_09_11_15_50_26_L2W.nc:rhos_826 
varname     : rhos_826 (Surface reflectance) 
name        : rhos_826 
unit        :        1 

HTH, I only see constant values in the tif or nc, though. (in this case GDAL has heuristics to determine the UTM zonal raster georeferencing in the netcdf, from the rhos_826#grid_mapping=transverse_mercator but I'm not sure how it gets the extent from that)

[mdsumner]

ah no, it's entirely standard, rectlinear x and y and grid_mapping:

        double transverse_mercator ;
                transverse_mercator:crs_wkt = "PROJCRS[\"WGS 84 / UTM zone 18N\",BASEGEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 18N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-75,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"easting\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"northing\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",32618]]" ;
                transverse_mercator:semi_major_axis = 6378137. ;
                transverse_mercator:semi_minor_axis = 6356752.31424518 ;
                transverse_mercator:inverse_flattening = 298.257223563 ;
                transverse_mercator:reference_ellipsoid_name = "WGS 84" ;
                transverse_mercator:longitude_of_prime_meridian = 0. ;
                transverse_mercator:prime_meridian_name = "Greenwich" ;
                transverse_mercator:geographic_crs_name = "WGS 84" ;
                transverse_mercator:horizontal_datum_name = "World Geodetic System 1984" ;
                transverse_mercator:projected_crs_name = "WGS 84 / UTM zone 18N" ;
                transverse_mercator:grid_mapping_name = "transverse_mercator" ;
                transverse_mercator:latitude_of_projection_origin = 0. ;
                transverse_mercator:longitude_of_central_meridian = -75. ;
                transverse_mercator:false_easting = 500000. ;
                transverse_mercator:false_northing = 0. ;
                transverse_mercator:scale_factor_at_central_meridian = 0.9996 ;
        double x(x) ;
                x:standard_name = "projection_x_coordinate" ;
                x:long_name = "x coordinate of projection" ;
                x:units = "m" ;
        double y(y) ;
                y:standard_name = "projection_y_coordinate" ;
                y:long_name = "y coordinate of projection" ;
                y:units = "m" ;

I may have missed something above, but certainly there's no GDAL issue here, with terra you may need to isolate to the NetCDF driver and pick the subdataset of interest (there are various options with DRIVER:file:sds or vrt://file?if=NetCDF&sd_name=sds for example).

[acolite]

I got this to work in a way that suprised me a bit, thanks for the example! (I can't speak on [@rhijmans] HTH, I only see constant values in the tif or nc, though. (in this case GDAL has heuristics to determine the UTM zonal raster georeferencing in the netcdf, from the rhos_826#grid_mapping=transverse_mercator but I'm not sure how it gets the extent from that)

The constant values are OK, as @arthurdegrandpre had to blank the image for licensing requirements. It does wonders for the compression though. 😄

[arthurdegrandpre]

Thanks Quiten, Michael, for your time, testing and inputs. I also do not have issues with the provided L2W products... I'm getting more and more convinced that it's just a hard to see user error at this point. Here is the complete output from my Acolite run (except the DEM), maybe you can see something that I am missing.

acolite_output.zip (It appears that I was mistaken, the issue appears between L1R and L2R, not from L2R to L2W, sorry about that)

Best, Arthur