0

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.

villesci
  • 5
  • 2
  • Not 100 % sure, but if I understood [tidyverse/ggplot2#3025](https://github.com/tidyverse/ggplot2/issues/3025#issuecomment-454373638) correctly, you should make use of `geom_raster()`. – dimfalk Aug 01 '22 at 22:55
  • Ah that is good to know! My cell resolution is not equal, so I think I should use geom_tile() for both df and raster plots. Using geom_tile() for the second plot actually helps a little (shifts closer to shapefile boundaries), but still not getting perfect line-up between figure1 and figure2 – villesci Aug 01 '22 at 23:41
  • Wait. Your cell resolution is equal - that's the whole point of using raster data. ERDDAP documentation says otherwise, I know: `geospatial_lat_resolution = 0.049999999999999996` and `geospatial_lon_resolution = 0.05000000000000001`, but I strongly assume this is some kind of imprecision of floating point data. Resolution is 0.05° in x and y direction. – dimfalk Aug 02 '22 at 22:29
  • Unfortunately your upper example as well as the one from `heatwaveR` vignette is failing, but you can have a glimpse in the dataset making use of the form provided - the result looks sth like [this](https://upwell.pfeg.noaa.gov/erddap/griddap/NOAA_DHW_Lon0360.htmlTable?CRW_SST%5B(2022-08-01T12:00:00Z):1:(2022-08-01T12:00:00Z)%5D%5B(3.00):1:(3.5)%5D%5B(0):1:(1)%5D). You can simply import the pasted table using `data <- readr::read_table()`, create a sf object via `sf::st_as_sf(data, coords = c("latitude", "longitude"), crs = sf::st_crs(4326))` and plot the result - a perfectly aligned grid. – dimfalk Aug 02 '22 at 22:35
  • 1
    Apologies for that - I did figure out that the reason for the shift in the plots was due to the extent in my raster - I calculated the minimum/maximum longitude and latitude from the SST data I obtained rather than using the extent I actually desired - the extracted values were 0.025 degrees off, which resulted in a shifted raster after reasterization. – villesci Aug 03 '22 at 16:39

0 Answers0