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.
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")