R stars: st_as_stars/read_stars kills x&y raster values; how to convert to stars?

Question

I have a raster read in & given a CRS with:

x <- raster::raster("originalobject")
dataCRS <- readRDS(paste0(crsloc, "CRS.Rds"))
raster::crs(x) <- dataCRS

It has an x extent of roughly -+29000, CRS is AEQD, 3857. Its data plot fine. I want to convert it to a stars object to plot it as a surface in a ggmap plot using geom_stars.

However if I try to convert the raster object to a stars object:

surfaceUD <- stars::st_as_stars(x)

x & y have wkt-style 'input' crs attributes but no unknown 'wkt' attributes, and x & y values are both NULL. The data in the raster ("All_Rasters_Scaled_Weighted_UDScaled") are still present:

[]

If I instead convert and set a CRS:

surfaceUD <- stars::st_as_stars(x) %>% sf::st_set_crs(3857)

x & y have 'input' crs attributes as "EPSG:3857", possibly having lost the lat/lon offsets? The 'wkt' attributes are now populated, and x & y values are both still NULL. The data in the raster are still present.

[]

Full x attribute info in case it's useful:

attr(surfaceUD, "dimensions")$x

$from 1

$to 293

$offset -29238.71

$delta 200

$refsys Coordinate Reference System: User input: EPSG:3857 wkt:

PROJCRS["WGS 84 / Pseudo-Mercator",

BASEGEOGCRS["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]],
    ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
    METHOD["Popular Visualisation Pseudo Mercator",
        ID["EPSG",1024]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["False easting",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]],
CS[Cartesian,2],
    AXIS["easting (X)",east,
        ORDER[3],
        LENGTHUNIT["metre",1]],
    AXIS["northing (Y)",north,
        ORDER[3],
        LENGTHUNIT["metre",1]],
USAGE[
    SCOPE["Web mapping and visualisation."],
    AREA["World between 85.06°S and 85.06°N."],
    BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]

$point NA

$values NULL

Possibly the $offset is the same as +lat_0=25.6871? I'm not sure that the offset is a problem though.

I also tried read_stars on "originalobject" but that had the same effect as surfaceUD <- stars::st_as_stars(x).

The stars intro says read_stars "will cause a loss of certain properties (cell size, reference system, vector geometries)". Is this the cause?

Is there any way to convert from raster to stars without losing all this vital info?

Thanks!

Answer 1

You can't use the CRS format from the sp package in st_set_crs from the sp package, you have to use crs from the raster package. raster::crs(x) <- dataCRS can accept a sp::CRS object, but sp::st_set_crs CAN'T accept a sp::CRS object.

Incredibly unintuitively, you have to use the proj4string minor function of raster::CRS:

x <- stars::read_stars(x) # x is an asc raster
dataCRS <- readRDS(paste0(crsloc, "CRS.Rds")) # load CRS from file
sf::st_crs(x) <- proj4string(dataCRS) # set CRS