problem with EPSG registry and GeoTIFF keys

[wiesehahn]

I try to use vapour to query a COG but I get a warning message and struggle to solve this. I guess its not related to vapour or gdal but more about the underlying data.

When I query the data

fileadress <- "https://dop20-rgb.opengeodata.lgln.niedersachsen.de/324965790/2019-04-07/dop20rgb_32_496_5790_2_ni_2019-04-07.tif"

cog.url <- paste0("/vsicurl/", fileadress)
ras <- terra::rast(cog.url)

roi <- sf::st_bbox(img1)
prj <- sf::st_crs(img1)$wkt
dim <-  c(500,500)

vals <- vapour::vapour_warp_raster(cog.url, extent = roi, dimension = dim, projection = prj)

I get the following warnings

proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
Warning 1: The definition of geographic CRS EPSG:4258 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

But I am not really getting the point here. Is this related to vapour or should I better switch to SO?

[wiesehahn]

BTW: gdalinfo for the image returns:

Size is 10000, 10000
Coordinate System is:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6┬░E and 12┬░E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
        BBOX[38.76,6,84.33,12.01]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
Origin = (552000.000000000000000,5726000.000000000000000)
Pixel Size = (0.200000000000000,-0.200000000000000)
Metadata:
  AREA_OR_POINT=Area
  OVR_RESAMPLING_ALG=BILINEAR
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  JPEGTABLESMODE=1
  JPEG_QUALITY=75
  LAYOUT=COG
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (  552000.000, 5726000.000) (  9d45' 7.76"E, 51d40'57.20"N)
Lower Left  (  552000.000, 5724000.000) (  9d45' 6.69"E, 51d39'52.47"N)
Upper Right (  554000.000, 5726000.000) (  9d46'51.90"E, 51d40'56.52"N)
Lower Right (  554000.000, 5724000.000) (  9d46'50.78"E, 51d39'51.79"N)
Center      (  553000.000, 5725000.000) (  9d45'59.28"E, 51d40'24.50"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313

[mdsumner]

Your code is incomplete, I can't reproduce the message with

u <- "https://dop20-rgb.opengeodata.lgln.niedersachsen.de/324965790/2019-04-07/dop20rgb_32_496_5790_2_ni_2019-04-07.tif"

cog.url <- file.path("/vsicurl", u)
info <- vapour::vapour_raster_info(cog.url)
roi <- info$extent
prj <- info$projection
dim <-  c(500,500)

vals <- vapour::vapour_warp_raster(cog.url, extent = roi, dimension = dim, projection = prj)

Need your sessionInfo() at least to explore what version/s system etc.

Also, sorry but extent and roi are incompatible - extent is xmin, xmax, ymin, ymax which you can get from sf with

roi <- sf::st_bbox(img1)[c(1, 3, 2, 4)]

fwiw, I can visualize the data like this

EDIT: orient the array correctly (this is the same as matrix( , dim[2], byrow = TRUE))

ximage::ximage(aperm(array(unlist(lapply(vals, as.integer)),  c(dim, 3)), c(2, 1, 3)), extent = roi, asp = 1)

or similarly with terra

img <- terra::rast(terra::ext(roi), ncols = dim[1], nrows = dim[2], vals = array(unlist(lapply(vals, as.integer)), c(dim, 3)), nlyrs = 3, crs = prj)
terra::plot(img)

I got some osmdata to check it lines up ok.

[wiesehahn]

Thanks for your quick reply. I got invalid data when running my code and thought it was related to the warning, but the actual reason was that the extent was invalid due to the different order in sf::st_bbox(). Applying your suggestion works!

So the wrning is unrelated to my initial problem, however it is still returned, e.g. by info <- vapour::vapour_raster_info(cog.url)

proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
Warning 1: The definition of geographic CRS EPSG:4258 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.2.2 tools_4.2.2    Rcpp_1.0.9     vapour_0.9.3   jsonlite_1.8.4

[mdsumner]

ah ok thank you, I see how to reproduce now too, I've not updated on windows recently 🙏

[mdsumner]

actually, no I don't know - it works fine here. If the data isn't obviously misaligned then I wouldn't worry about it.

u <- "https://dop20-rgb.opengeodata.lgln.niedersachsen.de/324965790/2019-04-07/dop20rgb_32_496_5790_2_ni_2019-04-07.tif"

 cog.url <- file.path("/vsicurl", u)
 info <- vapour::vapour_raster_info(cog.url)
 roi <- info$extent
 prj <- info$projection
 dim <-  c(500,500)
vals <- vapour::vapour_warp_raster(cog.url, extent = roi, dimension = dim, projection = prj)
> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
[3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.2.2 tools_4.2.2    parallel_4.2.2 Rcpp_1.0.9     vapour_0.9.3  
[6] jsonlite_1.8.4

I wonder, do you see the same message if you use othe package, say

terra::rast(cog.url) * 1

the *1 just ensures it pulls all the data locally, it shouldn't make a difference tho.

[wiesehahn]

I am getting only the first messages when reading with terra (doesnt matter if I pull the data locally or not)

> t1 <- terra::rast(cog.url) * 1
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
> t1 <- terra::rast(cog.url)
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found

but I get the entire warning when reading with stars

> t2 <- stars::read_stars(cog.url, proxy=FALSE)
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
Warning message:
In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
  GDAL Message 1: The definition of geographic CRS EPSG:4258 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

[mdsumner]

ok, at least it's not "my fault" alone :)

I don't have any capacity to chase this atm but I will if I get a chance. Happy for this to stay open until we can clear it. As I said, if it's working and not obviously misaligned spatially then I'd be comfortable ignore the warnings. If very precise measurement is crucial to you, or you notice a discrepancy that really matters in the data that's a different matter of course.