Projection issues after rasterize

73 views Asked by At

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)

Figure 1 output

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)

Figure2 output

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.

0

There are 0 answers