GDAL Message reading a tif file written with terra

[fgoerlich]

Hello,

A tif file written with terra is perfectly read (and manipulated) with this pkg:

terra::rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> class       : SpatRaster 
#> dimensions  : 1494, 1601, 1  (nrow, ncol, nlyr)
#> resolution  : 5, 5  (x, y)
#> extent      : 500497.5, 508502.5, 3901198, 3908668  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
#> source      : PNOA_MDT05_CodMuni_52001_buffer100.tif 
#> name        : PNOA_MDT05_ETRS89_HU30_1111_COR 
#> min value   :                          -4.534 
#> max value   :                         362.375 

but I get a GDAL warning if it is read with stars:

stars::read_stars("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>                                      Min. 1st Qu. Median     Mean 3rd Qu.    Max.
#> PNOA_MDT05_CodMuni_52001_buffe...  -0.797       0      0 75.20999 156.875 203.422
#> dimension(s):
#>   from   to  offset delta                refsys point values x/y
#> x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE   NULL [x]
#> y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE   NULL [y]
#> Warning messages:
#> 1: 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.
#> 2: 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.

I am not sure what this means exactly and how to proceed.
I am trying to move from terra to stars (mainly for compatibility with sf) so I have a lot of files written (and processed) with terra.

The (zipped) file is attached.

Thanks a lot!
PNOA_MDT05_CodMuni_52001_buffer100.zip

[edzer]

Doesn't happen here; can you pls report the output of terra::gdal() and sf::sf_extSoftVersion() as well as sessionInfo()?

[fgoerlich]

terra version 1.6.17

terra::gdal()
[1] "3.4.3"

sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.9.3"        "3.5.2"        "8.2.1"         "true"         "true"        "8.2.1" 

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

Matrix products: default

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

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

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         rstudioapi_0.14    magrittr_2.0.3     units_0.8-0        tidyselect_1.2.0   here_1.0.1         R6_2.5.1           rlang_1.0.6        fansi_1.0.3       
[10] dplyr_1.0.10       tools_4.2.2        tictoc_1.1         parallel_4.2.2     grid_4.2.2         utf8_1.2.2         KernSmooth_2.23-20 terra_1.6-17       cli_3.4.1         
[19] e1071_1.7-12       DBI_1.1.3          class_7.3-20       abind_1.4-5        assertthat_0.2.1   readxl_1.4.1       rprojroot_2.0.3    lwgeom_0.2-9       tibble_3.1.8      
[28] lifecycle_1.0.3    sf_1.0-9           stars_0.5-6        vctrs_0.5.1        codetools_0.2-18   glue_1.6.2         proxy_0.4-27       pillar_1.8.1       compiler_4.2.2    
[37] cellranger_1.1.0   generics_0.1.3     classInt_0.4-8     writexl_1.4.1      pkgconfig_2.0.3 

I have been unable to update the just released version of terra, 1.6-41:
installation of package ‘terra’ had non-zero exit status

[edzer]

This might be due to the differences in the respective GDAL versions used (I use 3.4.3 for both sf and terra); please try again when you've managed to install terra 1.6-41 - this should become available as binary package on CRAN within a few days. Cc: @rhijmans

[fgoerlich]

I have updated to terra version 1.6-41 which uses

terra::gdal()
[1] "3.5.2"

but things have worsened.
Reading the file from terra

terra::rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> class       : SpatRaster 
#> dimensions  : 1494, 1601, 1  (nrow, ncol, nlyr)
#> resolution  : 5, 5  (x, y)
#> extent      : 500497.5, 508502.5, 3901198, 3908668  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
#> source      : PNOA_MDT05_CodMuni_52001_buffer100.tif 
#> name        : PNOA_MDT05_ETRS89_HU30_1111_COR 
#> min value   :                          -4.534 
#> max value   :                         362.375 

while reading from stars I get the same output as before:

stars::read_stars("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>                                      Min. 1st Qu. Median     Mean 3rd Qu.    Max.
#> PNOA_MDT05_CodMuni_52001_buffe...  -0.797       0      0 75.20999 156.875 203.422
#> dimension(s):
#>   from   to  offset delta                refsys point x/y
#> x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE [x]
#> y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE [y]
#> Warning messages:
#> 1: 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.
#> 2: 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.

Should I open an issue in terra?

[edzer]

You can, for me it's not clear whether this is an issue with stars, terra, the .tif file, or some of the upstream libraries.

[fgoerlich]

You don´t see the issue with GDAL 3.4.3 in both terra and stars.
I see the issue with with GDAL 3.5.2 in both terra and stars, but I don´t see it with GDAL 3.4.3 in terra.
It seems an issue related with GDAL 3.5.2.
Is there anyway in which I can use GDAL 3.4.3 for both libraries.

[rhijmans]

It is a bit mysterious why you get these messages. Perhaps it is because you are using an old version of Windows (8.1).
I do not see them on Windows 10.

The message about an inconsistency in the CRS is odd. I see

library(terra)
r = rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
x = rast(crs="epsg:25830")
crs(r) == crs(x)
#[1] TRUE

Perhaps something went wrong with your installation. Or you are perhaps you are loading an old session?

[fgoerlich]

Yes, you are right. I am on Windows 8.1 on this machine. Is that the problem?
I am not loading an old session.
I may try the instalation again.

With your example I get:

> library(terra)
terra 1.6.41
> r = rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
> x = rast(crs="epsg:25830")
> crs(r) == crs(x)
[1] TRUE

Is it safe to use, or sholud I degrade terra on this machine?
It is difficult to undesrtand why the OS causes this with a more recent version of GDAL but not with the other.
Thanks a lot!

[rhijmans]

I do not know if this is related to the version of Windows. We also do not know for sure if this is related to the version of GDAL; but it could be. It looks like an incomplete installation; or a conflict between multiple installations of PROJ. But that is not something you would expect on Windows. Whether it matters I cannot say. Probably, but only you can find out. E.g. by trying to project the raster to another crs.

[fgoerlich]

I still see the problem in Windows 10 with GDAL 3.5.2, in both terra and stars

packageVersion("terra")
#> [1] ‘1.6.41’
terra::gdal()
#> [1] "3.5.2"
packageVersion("stars")
#> [1] ‘0.6.0’
sf::sf_extSoftVersion()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
#>        "3.9.3"        "3.5.2"        "8.2.1"         "true"         "true"        "8.2.1" 
> terra::rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> class       : SpatRaster 
#> dimensions  : 1494, 1601, 1  (nrow, ncol, nlyr)
#> resolution  : 5, 5  (x, y)
#> extent      : 500497.5, 508502.5, 3901198, 3908668  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
#> source      : PNOA_MDT05_CodMuni_52001_buffer100.tif 
#> name        : PNOA_MDT05_ETRS89_HU30_1111_COR 
#> min value   :                          -4.534 
#> max value   :                         362.375 
stars::read_stars("PNOA_MDT05_CodMuni_52001_buffer100.tif")
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> proj_create_from_database: datum not found
#> proj_create_from_database: ellipsoid not found
#> proj_create_from_database: prime meridian not found
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>                                      Min. 1st Qu. Median     Mean 3rd Qu.    Max.
#> PNOA_MDT05_CodMuni_52001_buffe...  -0.797       0      0 75.20999 156.875 203.422
#> dimension(s):
#>   from   to  offset delta                refsys point x/y
#> x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE [x]
#> y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE [y]
#> Warning messages:
#> 1: 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.
#> 2: 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.
sessionInfo()
#> R version 4.2.2 (2022-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#> 
#> Matrix products: default
#> 
#> locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8    LC_MONETARY=Spanish_Spain.utf8
[4] LC_NUMERIC=C                   LC_TIME=Spanish_Spain.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] reprex_2.0.2
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.9         highr_0.9          pillar_1.8.1       compiler_4.2.2     class_7.3-20      
#>  [6] tools_4.2.2        digest_0.6.30      evaluate_0.18      lifecycle_1.0.3    tibble_3.1.8      
#> [11] pkgconfig_2.0.3    rlang_1.0.6        DBI_1.1.3          cli_3.4.1          rstudioapi_0.14   
#> [16] yaml_2.3.6         writexl_1.4.1      parallel_4.2.2     xfun_0.34          fastmap_1.1.0     
#> [21] e1071_1.7-12       terra_1.6-41       withr_2.5.0        dplyr_1.0.10       knitr_1.41        
#> [26] generics_0.1.3     fs_1.5.2           vctrs_0.5.0        tictoc_1.1         classInt_0.4-8    
#> [31] grid_4.2.2         tidyselect_1.2.0   glue_1.6.2         sf_1.0-9           R6_2.5.1          
#> [36] processx_3.8.0     fansi_1.0.3        rmarkdown_2.18     callr_3.7.3        clipr_0.8.0       
#> [41] magrittr_2.0.3     ps_1.7.2           htmltools_0.5.3    codetools_0.2-18   stars_0.6-0       
#> [46] units_0.8-0        assertthat_0.2.1   abind_1.4-5        KernSmooth_2.23-20 utf8_1.2.2        
#> [51] proxy_0.4-27       lwgeom_0.2-9      

^{Created on 2022-11-23 with reprex v2.0.2}

[edzer]

You said that this tif file was written with terra. Does it also persist if you write it with updated terra, or with stars?

[fgoerlich]

Yes, I see the same issue if I write the file with the updated version of Terra.
Lets see the issue from first principles.
This is the original issue, reading an ESRI Ascii file and writte it with Terra 1.6-17:

> packageVersion("terra")
[1] ‘1.6.17’
> terra::gdal()
[1] "3.4.3"
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.9.3"        "3.5.2"        "8.2.1"         "true"         "true"        "8.2.1" 
> # Reading the original ascci-grid file and writting with terra
> r <- terra::rast("PNOA_MDT05_ETRS89_HU30_1111_COR.asc")
> terra::crs(r) <- "EPSG:25830"
> Melilla <- sf::read_sf("Melilla.gpkg")
> r <- terra::crop(r, terra::vect(sf::st_buffer(Melilla, dist = 100)))
> terra::writeRaster(r, "written_terra_1-6-17.tif")
> # Now I read the written file
> # terra 1.6.17 - GDAL 3.4.3 is ok
> terra::rast("written_terra_1-6-17.tif")
class       : SpatRaster 
dimensions  : 1494, 1601, 1  (nrow, ncol, nlyr)
resolution  : 5, 5  (x, y)
extent      : 500497.5, 508502.5, 3901198, 3908668  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
source      : written_terra_1-6-17.tif 
name        : PNOA_MDT05_ETRS89_HU30_1111_COR 
min value   :                          -4.534 
max value   :                         362.375 
> # stars GDAL 3.5.2 I get the warnings
> stars::read_stars("written_terra_1-6-17.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                            Min. 1st Qu. Median     Mean 3rd Qu.    Max.
written_terra_1-6-17.tif  -0.797       0      0 75.20999 156.875 203.422
dimension(s):
  from   to  offset delta                refsys point x/y
x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE [x]
y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE [y]
Warning messages:
1: 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.
2: 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.
> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8    LC_MONETARY=Spanish_Spain.utf8
[4] LC_NUMERIC=C                   LC_TIME=Spanish_Spain.utf8    

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

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         rstudioapi_0.14    magrittr_2.0.3     units_0.8-0        tidyselect_1.2.0   R6_2.5.1          
 [7] rlang_1.0.6        fansi_1.0.3        dplyr_1.0.10       tools_4.2.2        parallel_4.2.2     grid_4.2.2        
[13] KernSmooth_2.23-20 utf8_1.2.2         terra_1.6-17       cli_3.4.1          e1071_1.7-12       DBI_1.1.3         
[19] class_7.3-20       abind_1.4-5        assertthat_0.2.1   lwgeom_0.2-10      tibble_3.1.8       lifecycle_1.0.3   
[25] sf_1.0-9           stars_0.6-0        vctrs_0.5.1        codetools_0.2-18   glue_1.6.2         proxy_0.4-27      
[31] compiler_4.2.2     pillar_1.8.1       generics_0.1.3     classInt_0.4-8     pkgconfig_2.0.3

The tif file created is identical to the one attached at the begining of the discussion.

Now with Terra 1.6-41. Now I see the issue even when writing the file!

> packageVersion("terra")
[1] ‘1.6.41’
> terra::gdal()
[1] "3.5.2"
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.9.3"        "3.5.2"        "8.2.1"         "true"         "true"        "8.2.1" 
> # Reading the original ascci-grid file and writting with terra
> r <- terra::rast("PNOA_MDT05_ETRS89_HU30_1111_COR.asc")
> terra::crs(r) <- "EPSG:25830"
> Melilla <- sf::read_sf("Melilla.gpkg")
> r <- terra::crop(r, terra::vect(sf::st_buffer(Melilla, dist = 100)))
> terra::writeRaster(r, "written_terra_1-6-41.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
> # Now I read the written file
> # terra 1.6.41 - GDAL 3.4.3 is ok
> terra::rast("written_terra_1-6-41.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
class       : SpatRaster 
dimensions  : 1494, 1601, 1  (nrow, ncol, nlyr)
resolution  : 5, 5  (x, y)
extent      : 500497.5, 508502.5, 3901198, 3908668  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
source      : written_terra_1-6-41.tif 
name        : PNOA_MDT05_ETRS89_HU30_1111_COR 
min value   :                          -4.534 
max value   :                         362.375 
> # stars GDAL 3.5.2 I get the warnings
> stars::read_stars("written_terra_1-6-41.tif")
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                            Min. 1st Qu. Median     Mean 3rd Qu.    Max.
written_terra_1-6-41.tif  -0.797       0      0 75.20999 156.875 203.422
dimension(s):
  from   to  offset delta                refsys point x/y
x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE [x]
y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE [y]
Warning messages:
1: 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.
2: 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.
> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8    LC_MONETARY=Spanish_Spain.utf8
[4] LC_NUMERIC=C                   LC_TIME=Spanish_Spain.utf8    

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

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         rstudioapi_0.14    magrittr_2.0.3     units_0.8-0        tidyselect_1.2.0   R6_2.5.1          
 [7] rlang_1.0.6        fansi_1.0.3        dplyr_1.0.10       tools_4.2.2        parallel_4.2.2     grid_4.2.2        
[13] KernSmooth_2.23-20 utf8_1.2.2         terra_1.6-41       cli_3.4.1          e1071_1.7-12       DBI_1.1.3         
[19] class_7.3-20       abind_1.4-5        assertthat_0.2.1   lwgeom_0.2-10      tibble_3.1.8       lifecycle_1.0.3   
[25] sf_1.0-9           stars_0.6-0        vctrs_0.5.1        codetools_0.2-18   glue_1.6.2         proxy_0.4-27      
[31] compiler_4.2.2     pillar_1.8.1       generics_0.1.3     classInt_0.4-8     pkgconfig_2.0.3

Find attached the original file and the scripts.
issue_588.zip

[kadyb]

I think this problem is not related to {stars} / {terra} but to GDAL as I see the same messages directly in GDAL tools:

gdal_utils(util = "translate",
           source = "PNOA_MDT05_ETRS89_HU30_1111_COR.asc",
           destination = "test.tif",
           options = c("-a_srs", "EPSG:25830"))
gdal_utils(util = "info",
           source = "written_terra_1-6-17.tif")

Widnows 8.1; GDAL 3.5.2

[fgoerlich]

I see them also. Any chance to test it with GDAL 3.4.3?
I don´t see them with

terra::describe("written_terra_1-6-17.tif")

with

> packageVersion("terra")
[1] ‘1.6.17’
> terra::gdal()
[1] "3.4.3"

[edzer]

Am I right that this seems now more an upstream or a terra issue rather than a stars issue?

[kadyb]

Probably an upstream problem, because if we save this asc file with EPSG:25830 in {stars}, the same messages appear when loading.

Edit: @fgoerlich, maybe just reproject to EPSG:32630?

[fgoerlich]

sf::gdal_utils(util = "translate",
+   source = "PNOA_MDT05_ETRS89_HU30_1111_COR.asc",
+   destination = "test.tif",
+   options = c("-a_srs", "EPSG:32630"))

Saving the file with EPSG:32630 is ok. I don´t see any messages, either on writing or loading again.
Contrary to what happen when it is saved with EPSG:25830.

[kadyb]

In GDAL translate, you only overwritten the CRS. You should rather reproject raster in GDAL warp.

gdal_utils(util = "warp",
           source = "PNOA_MDT05_ETRS89_HU30_1111_COR.asc",
           destination = "test2.tif",
           options = c("-s_srs", "EPSG:25830",
                       "-t_srs", "EPSG:32630"))

[fgoerlich]

It is ok then. I don´t see any messages.
It is ok also if I do the oher way around.
Assuming the Ascii file is EPSG:32630 (which is not) and reproyecting to EPSG:25830, before saving.

sf::gdal_utils(util = "warp",
  source = "PNOA_MDT05_ETRS89_HU30_1111_COR.asc",
  destination = "test3.tif",
  options = c("-s_srs", "EPSG:32630",
              "-t_srs", "EPSG:25830"))

All this seems very odd.

[fgoerlich]

It is not necessary to save the file to see the messages.

> r <- terra::rast("PNOA_MDT05_CodMuni_52001_buffer100.tif")
> stars::st_as_stars(r)
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                                     Min. 1st Qu. Median     Mean 3rd Qu.    Max.
PNOA_MDT05_CodMuni_52001_buffe...  -0.797       0      0 75.20999 156.875 203.422
dimension(s):
  from   to  offset delta                refsys point x/y
x    1 1601  500498     5 ETRS89 / UTM zone 30N FALSE [x]
y    1 1494 3908668    -5 ETRS89 / UTM zone 30N FALSE [y]
Warning messages:
1: 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.
2: 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.

However If I project the raster (read with terra) to its own EPSG, then everything is ok:

> stars::st_as_stars(terra::project(r, "epsg:25830"))
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                                     Min. 1st Qu. Median     Mean 3rd Qu.    Max.
PNOA_MDT05_ETRS89_HU30_1111_CO...  -0.797       0      0 75.20999 156.875 203.422
dimension(s):
  from   to  offset delta                refsys x/y
x    1 1601  500498     5 ETRS89 / UTM zone 30N [x]
y    1 1494 3908668    -5 ETRS89 / UTM zone 30N [y]

Is this a workaround? or there is no problem at all?

[edzer]

I think there is no problem at all, since the file got read as ETRS89 / UTM zone 30N which is what was intended. If you add proxy = FALSE to read_stars you don't get the messages twice.

[fgoerlich]

Thanks! I still get the messages once with proxy = TRUE.
Is it safe to update to terra version 1.6-41, using GDAL 3.5.2?

[edzer]

Is it safe to update to terra version 1.6-41, using GDAL 3.5.2?

wrong place to ask, but in general you should always update to most recent CRAN versions.

[kadyb]

I think this is regression in GDAL or PROJ because today I encountered the same problem but for ESRI:102173 (ETRS 1989 Poland 1992). I'm sure it was working fine half year ago.

[kadyb]

@rouault, do you have any idea why these messages appear in stars and terra?

proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found

[rouault]

do you have any idea why these messages appear in stars and terra?

This slightly ring a bell of a bug I might have fixed some time ago, but can't find it. I don't reproduce this with recent GDAL versions when doing gdalinfo PNOA_MDT05_CodMuni_52001_buffer100.tif