Error cropping stars object to US census shape file

Question

I am working with a stars object containing temperature data from PRISM. The stars object contains daily temperature values for each latitude and longitude coordinate. However, I am trying to crop to a shape file of the contiguous U.S., and continue to encounter some issues.

Here is what I have attempted so far:

#Load packages
library(stars)
library(dplyr) 
library(raster)
library(terra)
library(tidyverse)

#Read in temperature file
tempdata <- readRDS("tempdata.rds")

tempdata has the following description/attributes:

stars object with 2 dimensions and 1 attribute
attribute(s):
                                  Min. 1st Qu. Median     Mean 3rd Qu.   Max.   NA's
PRISM_tmin_stable_4kmD2_202...  -9.449   9.476 16.621 15.18067  20.438 30.148 390874
dimension(s):
  from   to offset    delta refsys x/y
x    1 1405   -125  0.04167  NAD83 [x]
y    1  621  49.94 -0.04167  NAD83 [y]

#Confirm coordinate reference system- it is NAD83.
st_crs(tempdata)

#Read in shape file
shapefile <- st_read("shapefile.shp")

The shapefile has these attributes:

sing driver `ESRI Shapefile'
Simple feature collection with 33791 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -176.6967 ymin: -14.37378 xmax: 145.8305 ymax: 71.34132
Geodetic CRS:  NAD83

Then, I try to crop the larger prism file to the shape file:

#Crop to shape file
prism_cropped <- st_crop(tempdata, shapefile)

However, this returns the following error:

Warning message: In st_crop.stars(onedaytemp, shapefile) : st_crop: bounding boxes of x and y do not overlap

I do not know what to do at this point. There is definitely spatiotemporal data within the stars object, but for some reason it does not seem to be read properly/overlap with the shapefile, despite having the same coordinate reference system, and I am not sure how to fix it.

Answer 1

Looks like it is just a matter of converting the stars object to a spatRaster using the "terra" package. This code worked for me. Note, since the shapefile has over 30000 polygons, I just took the first five polygons from the shapefile so that it doesn't take forever to draw.

library(stars)
library(terra)
library(tmap)


#Read in temperature file
tempdata <- readRDS("PRISM_tmin_y20200701.rds") 
shapefile <- st_read("tl_2020_us_zcta520.shp")

temprast<-rast(tempdata)
croppedTemp<-crop(temprast, shapefile)

tmap_mode("plot")

tm_shape(croppedTemp)+
  tm_raster(palette="viridis")+
tm_shape(shapefile[1:5,])+
  tm_polygons(col=NA, border.col="red")