search

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
561 stars 94 forks source link

Remove proj_create_from_database warnings/messages #546

Closed paul-carteron closed 2 years ago

paul-carteron commented 2 years ago

Hi, When I'm using read_stars for a raster without coordinate there always warnings/messages and i'd like to remove them. I tried invisible(), suppressWarnings(), suppressMessages(), and the argument quiet = T, but nothing works, do you have any idea why ?

Here's a reprex :

library(stars)
#> Le chargement a nécessité le package : abind
#> Le chargement a nécessité le package : sf
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

url <- "https://wxs.ign.fr/altimetrie/geoportail/r/wms?version=1.3.0&request=GetMap&format=image/geotiff&layers=ELEVATION.ELEVATIONGRIDCOVERAGE&styles=&width=6&height=8&crs=EPSG:4326&bbox=47.79683,-4.375615,47.79859,-4.373898"

temp_file <- tempdir(check = TRUE)

download.file(url = url,
              method = "auto",
              mode = "wb",
              destfile = file.path(temp_file, "rast.tif"))

rast <- read_stars(file.path(temp_file, "rast.tif"), quiet = T)

# Messages I'd like to remove
#> 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

Created on 2022-07-10 by the reprex package (v2.0.1)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.0 (2022-04-22 ucrt) #> os Windows 10 x64 (build 19043) #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate French_France.utf8 #> ctype French_France.utf8 #> tz Europe/Paris #> date 2022-07-10 #> pandoc 2.17.1.1 @ C:/Program Files/RStudio/bin/quarto/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> abind * 1.4-5 2016-07-21 [1] CRAN (R 4.2.0) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.2.0) #> class 7.3-20 2022-01-16 [2] CRAN (R 4.2.0) #> classInt 0.4-7 2022-06-10 [1] CRAN (R 4.2.0) #> cli 3.3.0 2022-04-25 [1] CRAN (R 4.2.0) #> crayon 1.5.1 2022-03-26 [1] CRAN (R 4.2.0) #> DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.0) #> digest 0.6.29 2021-12-01 [1] CRAN (R 4.2.0) #> dplyr 1.0.9 2022-04-28 [1] CRAN (R 4.2.0) #> e1071 1.7-11 2022-06-07 [1] CRAN (R 4.2.0) #> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0) #> evaluate 0.15 2022-02-18 [1] CRAN (R 4.2.0) #> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0) #> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0) #> fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0) #> generics 0.1.2 2022-01-31 [1] CRAN (R 4.2.0) #> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0) #> highr 0.9 2021-04-16 [1] CRAN (R 4.2.0) #> htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.2.0) #> KernSmooth 2.23-20 2021-05-03 [2] CRAN (R 4.2.0) #> knitr 1.39 2022-04-26 [1] CRAN (R 4.2.0) #> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.2.0) #> lwgeom 0.2-8 2021-10-06 [1] CRAN (R 4.2.0) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0) #> pillar 1.7.0 2022-02-01 [1] CRAN (R 4.2.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0) #> proxy 0.4-27 2022-06-09 [1] CRAN (R 4.2.0) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.2.0) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0) #> Rcpp 1.0.8.3 2022-03-17 [1] CRAN (R 4.2.0) #> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.2.0) #> rlang 1.0.2 2022-03-04 [1] CRAN (R 4.2.0) #> rmarkdown 2.14 2022-04-25 [1] CRAN (R 4.2.0) #> rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.2.0) #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0) #> sf * 1.0-7 2022-03-07 [1] CRAN (R 4.2.0) #> stars * 0.5-5 2021-12-19 [1] CRAN (R 4.2.0) #> stringi 1.7.6 2021-11-29 [1] CRAN (R 4.2.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.2.0) #> tibble 3.1.7 2022-05-03 [1] CRAN (R 4.2.0) #> tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.2.0) #> units 0.8-0 2022-02-05 [1] CRAN (R 4.2.0) #> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0) #> vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.2.0) #> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0) #> xfun 0.31 2022-05-10 [1] CRAN (R 4.2.0) #> yaml 2.3.5 2022-02-21 [1] CRAN (R 4.2.0) #> #> [1] C:/Users/PAUL/AppData/Local/R/win-library/4.2 #> [2] C:/Program Files/R/R-4.2.0/library #> #> ────────────────────────────────────────────────────────────────────────────── ```
edzer commented 2 years ago

You can't silence them as they are emitted by the PROJ code, not by R. However they shouldn't be printed, it suggests there's something wrong with your PROJ or GDAL setup. What is your sessionInfo() after loading stars?

kadyb commented 2 years ago

Is the CRS of the input data correct?

library("stars")

tmp = tempfile(fileext = ".tif")
url = "https://wxs.ign.fr/altimetrie/geoportail/r/wms?version=1.3.0&request=GetMap&format=image/geotiff&layers=ELEVATION.ELEVATIONGRIDCOVERAGE&styles=&width=6&height=8&crs=EPSG:4326&bbox=47.79683,-4.375615,47.79859,-4.373898"
download.file(url, tmp, mode = "wb")
gdal_utils("info", tmp)
#> BOUNDCRS[
#>     SOURCECRS[
#>         GEOGCRS["unknown",
#>             DATUM["unnamed",
#>                 (...)

r = read_stars(tmp)
st_crs(r) = 4326
write_stars(r, tmp)
rr = read_stars(tmp) # no messages

In {terra} there are the same messages when loading.

paul-carteron commented 2 years ago

Hi @edzer, you're probably right because I test it an another PC and it works without showing warnings. Is there a way to clean GDAL and PROJ setup ? Here my sessionInfo() after loadding stars by the way :

Session info R version 4.2.0 (2022-04-22 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19043) Matrix products: default locale: [1] LC_COLLATE=French_France.utf8 LC_CTYPE=French_France.utf8 [3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C [5] LC_TIME=French_France.utf8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] stars_0.5-5 sf_1.0-7 abind_1.4-5 loaded via a namespace (and not attached): [1] Rcpp_1.0.9 rstudioapi_0.13 magrittr_2.0.3 [4] units_0.8-0 tidyselect_1.1.2 R6_2.5.1 [7] rlang_1.0.3 fansi_1.0.3 dplyr_1.0.9 [10] tools_4.2.0 parallel_4.2.0 grid_4.2.0 [13] KernSmooth_2.23-20 utf8_1.2.2 cli_3.3.0 [16] e1071_1.7-11 DBI_1.1.3 ellipsis_0.3.2 [19] class_7.3-20 assertthat_0.2.1 lwgeom_0.2-8 [22] tibble_3.1.7 lifecycle_1.0.1 crayon_1.5.1 [25] purrr_0.3.4 vctrs_0.4.1 glue_1.6.2 [28] proxy_0.4-27 compiler_4.2.0 pillar_1.7.0 [31] generics_0.1.3 classInt_0.4-7 pkgconfig_2.0.3

Hi @kadyb, you're right, the input data CRS are wrong but I can't help it. I'm doing the same thing as you inside my function. I'd just like to avoid the 6 lines of warnings because when I download several rasters, important informations displayed in the console are lost between warnings.

kadyb commented 2 years ago

Maybe as workaround you can overwrite the CRS of all files using gdal_translate and then load the data? ignore.stderr = TRUE suppresses the messages from PROJ. On Windows, you need to set an environment variable to run the GDAL from command line.

system(paste("gdal_translate", tmp, tmp2, "-a_srs EPSG:4326"),
       ignore.stderr = TRUE)

Edit: Or even better would be to use gdal_translate to download rasters from WMS and save them with corrected CRS (like here), but I haven't tested that.

paul-carteron commented 2 years ago

Hi @kabyb, thanks for the workaround but I before need to get familiar with GDAL.

For the moment, I try to understand from where the problem come from. I tried to reinstall stars from github, as source or binary. I tried same thing with sf. I tried to update GDAL and PROJ without sucess.

The only way I find to get around without a trick is by using older version of R 4.1.2 or R 4.1.0. which use GDAL 3.2.1 instead of GDAL 3.3.2 in last version of R.

edzer commented 2 years ago

If you use windows and install binary packages, updating GDAL or PROJ has no effect because they are part of the packages you install. If you install (compile) source packages (with rtools), this may have an effect (and solve the issue).

paul-carteron commented 2 years ago

I'm on windows and I compile with rtools42. I update PROJ with osgeo4w-setup, is this the right way to do it ?

At the end of source installation I can see those warnings if it can help : 

in the method for 'dbWriteTable' with the signature '"PostgreSQLConnection", "character", "sf"' : no definition of the class "PostgreSQLConnection" is found
in the method for 'dbDataType' with the signature '"PostgreSQLConnection", "sf"' : no definition of the class "PostgreSQLConnection" is found
in the method for 'coerce' with the signature '"Spatial", "sf"' : no definition of the class "Spatial" is found
in the method for 'coerce' with the signature '"Spatial", "sfc"' : no definition of the class "Spatial" is found
in the method for 'coerce' with the signature '"sf", "Spatial"' : no definition of the class "Spatial" is found
in the method for 'coerce' with the signature '"sfc", "Spatial"' : no definition of the class "Spatial" is found
in the method for 'coerce' with the signature '"XY", "Spatial"' : no definition of the class "Spatial" is found
in the method for 'coerce' with the signature '"crs", "CRS"' : no definition of the class "CRS" is found
in the method for 'coerce' with the signature '"sgbp", "sparseMatrix"' : no definition of the "sparseMatrix" class is found
edzer commented 2 years ago

The warnings you see are harmless.

I'm on windows and I compile with rtools42. I update PROJ with osgeo4w-setup, is this the right way to do it ?

I don't know much about windows; I would assume this is not the way to do this, because R 4.2 is linked against UCRT runtime, and I'd assume all libraries it links to should be as well. @rsbivand any idea?

rsbivand commented 2 years ago

If you mix binary libraries (RTools42 and OSGeo4W) built using different compilers and assumptions, random failures may appear. Crucially, RTools42 is up to speed on UCRT, but IIUC OSGeo4W is not.

Further, the messages indicate that the proj.db database if found is not the one used to install the package; This is not just about the PROJ library, but about the essential support files it uses. In GRASS GIS, this SQLite DB is asked for its version, and the version is checked against the install-time version.

rsbivand commented 2 years ago

The CRS being used is very unhelpful:

> st_crs(rast)
Coordinate Reference System:
  User input: unknown 
  wkt:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["unnamed",
                ELLIPSOID["unnamed",6378137,298.257223563,
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            ENSEMBLE["World Geodetic System 1984 ensemble",
                MEMBER["World Geodetic System 1984 (Transit)"],
                MEMBER["World Geodetic System 1984 (G730)"],
                MEMBER["World Geodetic System 1984 (G873)"],
                MEMBER["World Geodetic System 1984 (G1150)"],
                MEMBER["World Geodetic System 1984 (G1674)"],
                MEMBER["World Geodetic System 1984 (G1762)"],
                MEMBER["World Geodetic System 1984 (G2139)"],
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ENSEMBLEACCURACY[2.0]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            USAGE[
                SCOPE["Horizontal component of 3D system."],
                AREA["World."],
                BBOX[-90,-180,90,180]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

The source appears to be adding a +towgs84= key. So the data source may also be part of the problem.

paul-carteron commented 2 years ago

Hi @rsbivand, I don't get same crs as yours when I tried with R 4.2 as you can see in the reprex below. But, when I use R 4.1, I get same crs and no warnings are shown. I have the impression that the problem comes from the version of GDAL used (3.2.1 works but 3.4.3 doesn't).

library(stars)
#> Le chargement a nécessité le package : abind
#> Le chargement a nécessité le package : sf
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

url <- "https://wxs.ign.fr/altimetrie/geoportail/r/wms?version=1.3.0&request=GetMap&format=image/geotiff&layers=ELEVATION.ELEVATIONGRIDCOVERAGE&styles=&width=6&height=8&crs=EPSG:4326&bbox=47.79683,-4.375615,47.79859,-4.373898"

temp_file <- tempdir(check = TRUE)

download.file(url = url,
              method = "auto",
              mode = "wb",
              destfile = file.path(temp_file, "rast.tif"))

rast <- read_stars(file.path(temp_file, "rast.tif"), quiet = T)
st_crs(rast)
#> Coordinate Reference System:
#>   User input: unknown 
#>   wkt:
#> GEOGCRS["unknown",
#>     DATUM["unnamed",
#>         ELLIPSOID["unnamed",6378137,298.257223563,
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433,
#>             ID["EPSG",9122]]],
#>     CS[ellipsoidal,2],
#>         AXIS["latitude",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["longitude",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

Created on 2022-07-14 by the reprex package (v2.0.1)

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.1 (2022-06-23 ucrt) #> os Windows 10 x64 (build 19043) #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate French_France.utf8 #> ctype French_France.utf8 #> tz Europe/Paris #> date 2022-07-14 #> pandoc 2.17.1.1 @ C:/Program Files/RStudio/bin/quarto/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> abind * 1.4-5 2016-07-21 [1] CRAN (R 4.2.0) #> class 7.3-20 2022-01-16 [2] CRAN (R 4.2.1) #> classInt 0.4-7 2022-06-10 [1] CRAN (R 4.2.1) #> cli 3.3.0 2022-04-25 [1] CRAN (R 4.2.1) #> crayon 1.5.1 2022-03-26 [1] CRAN (R 4.2.1) #> DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.1) #> digest 0.6.29 2021-12-01 [1] CRAN (R 4.2.1) #> dplyr 1.0.9 2022-04-28 [1] CRAN (R 4.2.1) #> e1071 1.7-11 2022-06-07 [1] CRAN (R 4.2.1) #> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.1) #> evaluate 0.15 2022-02-18 [1] CRAN (R 4.2.1) #> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.1) #> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.1) #> fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.1) #> generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1) #> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.1) #> highr 0.9 2021-04-16 [1] CRAN (R 4.2.1) #> htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.2.1) #> KernSmooth 2.23-20 2021-05-03 [2] CRAN (R 4.2.1) #> knitr 1.39 2022-04-26 [1] CRAN (R 4.2.1) #> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.2.1) #> lwgeom 0.2-8 2021-10-06 [1] CRAN (R 4.2.1) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.1) #> pillar 1.7.0 2022-02-01 [1] CRAN (R 4.2.1) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.1) #> proxy 0.4-27 2022-06-09 [1] CRAN (R 4.2.1) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.2.1) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.1) #> Rcpp 1.0.9 2022-07-08 [1] CRAN (R 4.2.1) #> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.2.1) #> rlang 1.0.4 2022-07-12 [1] CRAN (R 4.2.1) #> rmarkdown 2.14 2022-04-25 [1] CRAN (R 4.2.1) #> rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.2.1) #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.1) #> sf * 1.0-8 2022-07-14 [1] CRAN (R 4.2.1) #> stars * 0.5-5 2021-12-19 [1] CRAN (R 4.2.1) #> stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.1) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.2.1) #> tibble 3.1.7 2022-05-03 [1] CRAN (R 4.2.1) #> tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.2.1) #> units 0.8-0 2022-02-05 [1] CRAN (R 4.2.1) #> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.1) #> vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.2.1) #> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.1) #> xfun 0.31 2022-05-10 [1] CRAN (R 4.2.1) #> yaml 2.3.5 2022-02-21 [1] CRAN (R 4.2.0) #> #> [1] C:/Users/PAUL/AppData/Local/R/win-library/4.2 #> [2] C:/Program Files/R/R-4.2.1/library #> #> ────────────────────────────────────────────────────────────────────────────── ```
kadyb commented 2 years ago

I have these setups and messages appear on both of them:

Windows 8.1; R 4.2: Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
Pop!_OS 20.04 LTS; R 4.2: Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

I think the easiest solution would be if you just overwrite the CRS of the rasters.

BTW: Isn't it strange that {stars} returns 6 messages, and {terra} and gdal_utils only return 3? Is there one duplicated hook from {stars} to PROJ?

edzer commented 2 years ago

stars might open the file twice if proxy is not given in read_stars().

rsbivand commented 2 years ago

My PROJ is 9.0.1, so gives datum ensembles introduced in 8. Nothing to do with GDAL. If newer PROJ and/or GDAL give different output, it is either because of fixed bugs, or newer versions of proj.db. The input file has a malformed CRS description.

edzer commented 2 years ago

The input file has a malformed CRS description.

I think that closes the issue?

paul-carteron commented 2 years ago

Yes, thank you all for your help and your time !