I am having issues when plotting SST data after rasterizing the original data file. I eventually want to extract SST from a series of points for each day over a period of time, and so I noticed a shift in the raster relative to coastline after I recieved a bunch of NAs after extraction. I suspect I messed up my projection/transformation, but everything seems to work.
Query to get 1 month of SST data as a tibble. From heatwaveR documentation. CRS for this data is EPSG:4326 according to ERDDAP. Resolution needed for later also from this step.
library(rerddap)
library(ggplot2)
library(dplyr)
library(raster)
library(rasterVis)
library(viridis)
library(purrr)
library(sf)
SST_sub_dl <- function(time_df){
SST_data_tas <- griddap(x = "NOAA_DHW_Lon0360",
url = "https://coastwatch.pfeg.noaa.gov/erddap/",
time = c(time_df$start, time_df$end),
latitude = c(-44.5, -38.5),
longitude = c(142, 150),
store=disk(),
fields = "CRW_SST")$data%>%
mutate(time = as.Date(stringr::str_remove(time, "T00:00:00Z"))) %>%
dplyr::rename(t = time, temp = CRW_SST) %>%
select(lon, lat, t, temp) %>%
na.omit()
}
#Assign years desired
tas_time <- data.frame(date_index = 1:2,
start = as.Date(c("2011-01-01")),
end = as.Date(c("2011-01-31")))
# The time this takes will vary greatly based on connection speed
#takes me 26 seconds
system.time(
SST_data_tas <- tas_time %>%
group_by(date_index) %>%
group_modify(~SST_sub_dl(.x)) %>%
ungroup() %>%
select(lon, lat, t, temp))
Now get Australia shapefile
aus.gadm<-getData("GADM",country="AUS",level=1, path=tempdir())
tas.gadm<-filter(aus.gadm,NAME_1=="Tasmania")
#clip function
gClip <- function(shp, bb){
if(class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))), "SpatialPolygons")
else b_poly <- as(extent(bb), "SpatialPolygons")
gIntersection(shp, b_poly, byid = T)
}
#clip to remove outer islands
b<-st_bbox(c(xmin=142,xmax=150,ymin=-44.5,ymax=-38.5),crs=4326)
tasss.gadm<-gClip(tas.gadm,b)
tasss.gadm.fort<-sf::st_as_sf(tasss.gadm,"+proj=longlat +datum=WGS84 +nodefs")
Here, I plot the SST data from the dataframe and get SST data and Tasmania shapefile lined up well.
figure1<-SST_data_tas %>%
filter(t == "2011-01-01") %>%
ggplot(aes(x = lon, y = lat)) +
geom_tile(aes(fill = temp)) +
scale_fill_viridis_c(na.value="transparent") +
coord_quickmap(expand = F) +
labs(x = NULL, y = NULL, fill = "SST (°C)") +
theme(legend.position = "bottom")+
geom_sf(data=tasss.gadm.fort,inherit.aes=F,fill=NA)
Of course, since I eventually want to extract SST values underneath points, I need to convert this to a raster stack, one layer for each day.
#create blank raster
r_tas_obj<- raster(xmn=min(SST_data_tas$lon), xmx=max(SST_data_tas$lon), ymn=min(SST_data_tas$lat),
ymx=max(SST_data_tas$lat),res=c(0.04999,0.05),crs=sf::st_crs(4326)[[2]])
#rasterize SST data into a rasterstack
system.time(SST_data_tas_stack <- SST_data_tas %>%
group_split(t) %>%
purrr::map(~rasterize(x=.x[,c("lon", "lat")],y=r_tas_obj,field=.x[,4],fun=mean))%>%
stack())
#plot
figure2<-SST_data_tas_stack$layer.1%>%
gplot()+
geom_raster(aes(x = x, y = y,fill=value),stat='identity') +
scale_fill_viridis_c(na.value="transparent") +
coord_quickmap(expand = F) +
labs(x = NULL, y = NULL, fill = "SST (°C)") +
theme(legend.position = "bottom")+
geom_sf(data=tasss.gadm.fort,inherit.aes=F,fill=NA)
#CRS match...
compareCRS(SST_data_tas_stack,tasss.gadm.fort)
What you will notice is that my raster result is slightly shifted compared to the original plotting of the data using geom_tile. I suspect this is due to an error at the rasterization process, but the CRS between my tasmania shapefile and SST raster match up.