Possible problem with TMAP and raster 3.3-7

[walts]

I had been successfully using the following R code until I updated to raster release 3.3-7 a few days ago

   library(raster)
   library(tmap)

   download.file(url = "https://drive.google.com/uc?export=download&id=1-C48DzgpuXwKwCzeeI57vrkel66YEntF",
              destfile = "houston_wards.jpg", mode="wb")

   ward_map <- raster::raster(x = "houston_wards.jpg",
                           crs= "+proj=longlat +datum=WGS84")  

   xmin(ward_map) <- -95.42946
   xmax(ward_map) <- -95.27887
   ymin(ward_map) <-  29.706078
   ymax(ward_map) <-  29.812857

   ward_map1 <- tm_shape(ward_map) + tm_raster(alpha = .5, palette = "Greys", legend.show = FALSE)
   print(ward_map1)

With TMAP 3.0, i got the following error:

"Error: Current projection of shape ward_map unknown and cannot be determined."

Using the development 3.1 version of TMAP, I now get the error

"stars_proxy object shown at 1109 by 901 cells.
Error: Current projection of shape ward_map unknown and cannot be determined."

Any thoughts?

[mtennekes]

This is probably an issue with CRS to crs conversion:

library(raster)
#> Loading required package: sp
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(stars)
#> Loading required package: abind

download.file(url = "https://drive.google.com/uc?export=download&id=1-C48DzgpuXwKwCzeeI57vrkel66YEntF",
              destfile = "houston_wards.jpg", mode="wb")

ward_map <- raster::raster(x = "houston_wards.jpg",
                           crs= "+proj=longlat +datum=WGS84")  

xmin(ward_map) <- -95.42946
xmax(ward_map) <- -95.27887
ymin(ward_map) <-  29.706078
ymax(ward_map) <-  29.812857

ward_stars = st_as_stars(ward_map)

st_crs(ward_stars)
#> Coordinate Reference System: NA

ward_stars2 = ward_stars %>% 
    st_set_crs("+proj=longlat +datum=WGS84")

st_crs(ward_stars2)
#> Coordinate Reference System:
#>   User input: +proj=longlat +datum=WGS84 
#>   wkt:
#> GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

The only difference is that user input is not transferred when converting a CRS to a crs:

st_crs(raster::crs("+proj=longlat +datum=WGS84"))
#> Coordinate Reference System:
#>   User input: unknown 
#>   wkt:
#> GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

st_crs("+proj=longlat +datum=WGS84")
#> Coordinate Reference System:
#>   User input: +proj=longlat +datum=WGS84 
#>   wkt:
#> GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

Any ideas @edzer?

[edzer]

Does tmap use the user input then?

[mtennekes]

No, but is does use st_as_stars to convert raster objects.
The problem seems to be

ward_stars = st_as_stars(ward_map)
st_crs(ward_stars)
#> Coordinate Reference System: NA

while ward_map has a valid CRS.

[edzer]

Ah, thanks - yes this is the (newer) use case where Raster objects are on-disk; stars now passes on raster's CRS, even if this is a nonsense CRS as in this case.

[mtennekes]

Thx! The CRS part has been fixed. I found another issue, namely that the extent is not converted as well.
I've just made a pull-request for that (r-spatial/stars#303).

With these updates to stars, the result is:

tm_shape(ward_map) + tm_rgb()