[gdal-dev] NetCDF and ESA Probav issue with corner coordimates

Sun Nov 3 10:00:56 PST 2019

GMT pays particular attention to this matter and detects grid vs pixel registration by checking the X,Y vector coordinates versus array sizes (for netCDF). It normally has no problem detecting the two cases.

Joaquim

From: gdal-dev <gdal-dev-bounces at lists.osgeo.org> On Behalf Of Ivan Lucena
Sent: Friday, November 1, 2019 1:18 PM
To: Andrew C Aitchison <andrew at aitchison.me.uk>; Richard Duivenvoorde <rdmailings at duif.net>
Cc: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] NetCDF and ESA Probav issue with corner coordimates

Thanks Andrew and Richard,

Very good examples.

I have struggled with that so many times when developing GDAL drivers or when dealing with data in formats like HDF4, HDF5, NetCDF.

Some formats have a pretty good easy way of telling which way to go about upperleft or center pixel coordinate. But that is not the case with "flexible" formats like the ones above. In some cases developer have to identify the data product in a hardcode way and apply what the documentation says about the data, not only in regards to pixel coordinate but other issues.

What is interesting in that case is that one can say that according to the documentation the driver is correct. The -180,80 coordinate is in the center of the pixel. But when putting that date to real use, that seems to be incorrect.

And that could be dangerous or at least inconvenient. In my case I have a Python code that does zonal statistic and when a geometry is converted to a little raster binary mask, it does not match the same area in the ProbaV data and the results is a crash in Fiji because ReadRaster goes overboard or lost of NDVI data in other regions. The half pixel shift grows as you go east.

In other case, a GIS user might find impossible to run a raster tool, because it will complain that all the raster sources that they are trying to combine does not match the same area. There are still some tool like that in the marked, unfortunately.

If I have some time, I will try to run the NetCDF driver in step-by-step debugging mode to see exactly how the crs:GeoTransform is (mis)handled.

The goal is to see if the data is wrong or if the driver is wrong (or if I am wrong) and to decided how to go from here. We might need to add some hardcoded decision into the NetCDF driver, which I would prefer not to do.

Best regards,

Ivan

________________________________
From: Andrew C Aitchison <andrew at aitchison.me.uk<mailto:andrew at aitchison.me.uk>>
Sent: Friday, November 1, 2019 7:24 AM
To: Richard Duivenvoorde <rdmailings at duif.net<mailto:rdmailings at duif.net>>
Cc: Ivan Lucena <ivan.lucena at outlook.com<mailto:ivan.lucena at outlook.com>>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org> <gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>>
Subject: Re: [gdal-dev] NetCDF and ESA Probav issue with corner coordimates

Does
   https://gis.stackexchange.com/questions/122670/is-there-a-standard-for-the-coordinates-of-pixels-in-georeferenced-rasters#122687
help you to understand the issue ?

On Fri, 1 Nov 2019, Richard Duivenvoorde wrote:

> Hi Ivan,
>
> I think I reported an issue like this some time ago in the QGIS issue
> tracker:
>
> https://issues.qgis.org/issues/20730
>
> I think(!), but I hope some netcdf guru jumps in, that the crux is how
> you define/interpret your data: is the value valid in the center of the
> grid (like a raster) or is it on the vertices (more like a mesh).
> It is all about interpretation of a datamodel: is it a grid or a raster.
>
> If I understand correctly in a mesh (and netcdf) is much more possible
> then in a raster: triangular meshes, non-regular etc etc:
>
> https://docs.qgis.org/3.4/en/docs/user_manual/working_with_mesh/mesh_properties.html
>
> I'm not sure who should be responsible for 'correcting/transforming'
> this: a netcdf which is written as if the values are on (for example)
> the upper-left corner vertices, while you actually mean the value is in
> the center of a (grid)cell... For a simple rectangular mesh, this seems
> feasible on the client side (see panoply example in first link), but for
> non-regular meshes/netcdf's this seems impossible client side (well
> actually that is not a raster/grid anymore)..
>
> In your case I would think that the data is actually mend to be a grid...
>
> Not sure if this is helpfull, as said I hope some netcdf guru jumps in.
> Maybe I'm oversimplifying here...
>
> Regards,
>
> Richard Duivenvoorde
>
> On 31/10/2019 22.06, Ivan Lucena wrote:
>>
>> Hi Folks,
>>
>> I Downloaded NetCDF files
>> from?http://land.copernicus.eu/global/products/ndvi and I am having
>> problems with the GeoTransform that the driver is getting.
>>
>> The image extents, as you can see in following gdalinfo report, goes
>> beyond the -180 and +180 longitude degrees (half resolution, half pixel)
>> but that's OK according to page 27 in the PDF documentation:
>> PROBAV-Products_User_Manual_v1.3.pdf
>> <https://earth.esa.int/documents/700255/1929094/PROBAV-Products_User_Manual_v1.3.pdf/fd5e30f4-5305-4d41-86aa-fab7b9568251>?
>>
>> But if I open that image on QGIS and overlay a vector (GADM from
>> gadm.org) or any other global layer it is clear that -180x80 should be
>> in the upper left corner of the image, not In the center. There is a
>> half pixel missing on the left and bottom of the image. This in in Fiji:
>>
>>
>>
>> The gdalinfo report shows that the crs:GeoTransform is?-180.0000000000
>> 0.0029761905 0.0 80.0000000000 0.0 -0.0029761905 but it seems like the
>> driver decided to take that to as center of the pixel based. Looking at
>> the NetCDF source I can see that there this possibility exists.
>>
>> Of course, I am not waiting for a solution. I just used gdal_translate
>> "- a_ullr" and fix that into the target format. It takes centuries to
>> converted that image to geotiff. But I am concerned that others might
>> get in trouble because of that.
>>
>> I tried to contact the data provider but did not received an answer yet.
>>
>> Does anybody has experience with that dataset and/or NetCDF that can
>> help figure out what is going on. It is a bug in the driver or an issue
>> with the data.
>>
>> Thanks,
>>
>> -
>> Ivan
>>
>> Here is the gdalinfo report:
>>
>> Don't worry about the warnings in regards to the pixel values. It is fine.
>>
>> $ gdalinfo netcdf:c_gls_NDVI300_201910110000_GLOBE_PROBAV_V1.0.1.nc
>> --debug on
>> GDAL_netCDF: driver detected file type=5, libnetcdf detected type=3
>> Warning 1: NetCDF driver detected file type=5, but libnetcdf detected
>> type=3
>> GDAL_netCDF: setting file type to 3, was 5
>> GDAL_netCDF: dim_count = 2
>> GDAL_netCDF: var_count = 4
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute _FillValue
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute missing_value
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute flag_values
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute valid_range
>> GDAL_netCDF:
>> =====
>> SetProjectionFromVar( 3)
>> *GDAL_netCDF: got grid_mapping crs*
>> GDAL_netCDF: bIsGdalFile=1 bIsGdalCfFile=0 bBottomUp=0
>> GDAL_netCDF: got spheroid from CF: (6378137.000000 , 298.257224)
>> GDAL_netCDF: set bBottomUp = 0 from Y axis
>> GDAL_netCDF: setting WKT from CF
>> GDAL_netCDF: SetProjection, WKT =
>> GEOGCS["unknown",DATUM["unknown",SPHEROID["Spheroid",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]]
>> 
>> GDAL_netCDF: xdim: 120960 nSpacingBegin: 3 nSpacingMiddle: 3
>> nSpacingLast: 3
>> GDAL_netCDF: ydim: 47040 nSpacingBegin: -3 nSpacingMiddle: -3
>> nSpacingLast: -3
>> GDAL_netCDF: setting WKT from GDAL
>> GDAL_netCDF: SetProjection, WKT = GEOGCS["WGS
>> 84",DATUM["WGS_1984",SPHEROID["WGS
>> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]
>> 
>> *GDAL_netCDF:
>> SetGeoTransform(-180.001488,0.002976,0.000000,80.001488,0.000000,-0.002976)
>> *
>> GDAL_netCDF: bGotGeogCS=1 bGotCfSRS=1 bGotCfGT=1 bGotGdalSRS=1 bGotGdalGT=0
>> Warning 1: Unsupported netCDF datatype (7), treat as Float32.
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute _FillValue
>> GDAL_netCDF: netcdf type=5 gdal type=6 signedByte=1
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute _FillValue
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute missing_value
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute flag_values
>> GDAL_netCDF: NCDFGetAttr unsupported type 7 for attribute valid_range
>> GDAL_netCDF: got add_offset=-0.07999999821186066, status=0
>> GDAL_netCDF: got scale_factor=0.004000000189989805, status=0
>> GDAL: GDALOpen(netcdf:c_gls_NDVI300_201910110000_GLOBE_PROBAV_V1.0.1.nc,
>> this=00000219E2D291F0) succeeds as netCDF.
>> Driver: netCDF/Network Common Data Format
>> GDAL: GDALDefaultOverviews::OverviewScan()
>> Files: none associated
>> Size is 120960, 47040
>> Coordinate System is:
>> GEOGCS["WGS 84",
>> ? ? DATUM["WGS_1984",
>> ? ? ? ? SPHEROID["WGS 84",6378137,298.257223563,
>> ? ? ? ? ? ? AUTHORITY["EPSG","7030"]],
>> ? ? ? ? TOWGS84[0,0,0,0,0,0,0],
>> ? ? ? ? AUTHORITY["EPSG","6326"]],
>> ? ? PRIMEM["Greenwich",0,
>> ? ? ? ? AUTHORITY["EPSG","8901"]],
>> ? ? UNIT["degree",0.0174532925199433,
>> ? ? ? ? AUTHORITY["EPSG","9108"]],
>> ? ? AUTHORITY["EPSG","4326"]]
>> *Origin = (-180.001488095238102,80.001488095238088)*
>> Pixel Size = (0.002976190476204,-0.002976190476190)
>> Metadata:
>> *? crs#GeoTransform=-180.0000000000 0.0029761905 0.0 80.0000000000 0.0
>> -0.0029761905*
>> ? crs#grid_mapping_name=latitude_longitude
>> ? crs#inverse_flattening=298.257223563
>> ? crs#longitude_of_prime_meridian=0
>> ? crs#long_name=coordinate reference system
>> ? crs#semi_major_axis=6378137
>> ? crs#spatial_ref=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
>> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]
>> 
>> ? crs#_CoordinateAxisTypes=GeoX GeoY
>> ? crs#_CoordinateTransformType=Projection
>> ? lat#axis=Y
>> ? lat#DIMENSION_LABELS=lat
>> ? lat#long_name=latitude
>> ? lat#standard_name=latitude
>> ? lat#units=degrees_north
>> ? lat#_CoordinateAxisType=Lat
>> ? lon#axis=X
>> ? lon#DIMENSION_LABELS=lon
>> ? lon#long_name=longitude
>> ? lon#standard_name=longitude
>> ? lon#units=degrees_east
>> ? lon#_CoordinateAxisType=Lon
>> ? NC_GLOBAL#archive_facility=VITO
>> ? NC_GLOBAL#Conventions=CF-1.6
>> ? NC_GLOBAL#copyright=Copernicus Service information 2019
>> ? NC_GLOBAL#history=Processing line NDVI: 2019-10-22
>> ?
>> NC_GLOBAL#identifier=urn:cgls:global:ndvi300_v1_333m:NDVI300_201910110000_GLOBE_PROBAV_V1.0.1
>> 
>> ? NC_GLOBAL#institution=VITO NV
>> ? NC_GLOBAL#long_name=Normalized Difference Vegetation Index
>> ? NC_GLOBAL#orbit_type=LEO
>> ? NC_GLOBAL#parent_identifier=urn:cgls:global:ndvi300_v1_333m
>> ? NC_GLOBAL#platform=Proba-V
>> ? NC_GLOBAL#processing_level=L3
>> ? NC_GLOBAL#processing_mode=Near Real Time
>> ? NC_GLOBAL#product_version=V1.0.1
>> ? NC_GLOBAL#references=http://land.copernicus.eu/global/products/ndvi
>> ? NC_GLOBAL#sensor=VEGETATION
>> ? NC_GLOBAL#source=Derived from EO satellite imagery
>> ? NC_GLOBAL#time_coverage_end=2019-10-20T23:59:59Z
>> ? NC_GLOBAL#time_coverage_start=2019-10-11T00:00:00Z
>> ? NC_GLOBAL#title=10-daily Normalized Difference Vegetation Index 333M:
>> GLOBE 2019-10-11T00:00:00Z
>> ? NDVI#add_offset=-0.079999998
>> ? NDVI#flag_meanings=Missing cloud snow sea background
>> ? NDVI#flag_values={}
>> ? NDVI#grid_mapping=crs
>> ? NDVI#long_name=Normalized Difference Vegetation Index 333M
>> ? NDVI#missing_value=
>> ? NDVI#scale_factor=0.0040000002
>> ? NDVI#standard_name=normalized_difference_vegetation_index
>> ? NDVI#units=
>> ? NDVI#valid_range={}
>> ? NDVI#_FillValue=
>> OGRCT: PROJ >= 4.8.0 features enabled
>> OGRCT: Using locale-safe proj version
>> OGRCT: Source: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
>> OGRCT: Target: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
>> Corner Coordinates:
>> *Upper Left ?(-180.0014881, ?80.0014881) (180d 0' 5.36"W, 80d 0' 5.36"N)*
>> *Lower Left ?(-180.0014881, -59.9985119) (180d 0' 5.36"W, 59d59'54.64"S)*
>> *Upper Right ( 179.9985119, ?80.0014881) (179d59'54.64"E, 80d 0' 5.36"N)*
>> *Lower Right ( 179.9985119, -59.9985119) (179d59'54.64"E, 59d59'54.64"S)*
>> Center ? ? ?( ?-0.0014881, ?10.0014881) ( ?0d 0' 5.36"W, 10d 0' 5.36"N)
>> Band 1 Block=120960x1 Type=Float32, ColorInterp=Undefined
>> ? NoData Value=0
>> ? Offset: -0.0799999982118607, ? Scale:0.00400000018998981
>> ? Metadata:
>> ? ? add_offset=-0.079999998
>> ? ? flag_meanings=Missing cloud snow sea background
>> ? ? flag_values={}
>> ? ? grid_mapping=crs
>> ? ? long_name=Normalized Difference Vegetation Index 333M
>> ? ? missing_value=
>> ? ? NETCDF_VARNAME=NDVI
>> ? ? scale_factor=0.0040000002
>> ? ? standard_name=normalized_difference_vegetation_index
>> ? ? units=
>> ? ? valid_range={}
>> ? ? _FillValue=
>> GDAL:
>> GDALClose(netcdf:c_gls_NDVI300_201910110000_GLOBE_PROBAV_V1.0.1.nc,
>> this=00000219E2D291F0)
>>
>>
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20191103/c1b6b0e4/attachment-0001.html>