Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ggmap rasters #20

Open
quantifish opened this issue Aug 5, 2023 · 5 comments
Open

ggmap rasters #20

quantifish opened this issue Aug 5, 2023 · 5 comments

Comments

@quantifish
Copy link
Owner

quantifish commented Aug 5, 2023

# https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-raster
library(ggplot2)
library(ggmap)
library(sf)
nc_map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google")
st_crs(nc_map)
# Coordinate Reference System: NA
# assume the coordinate refence system is 3857
plot(st_transform(nc_shp, crs = 3857)[1], bgMap = nc_map)
@quantifish
Copy link
Owner Author

register_google(key = key)

ggmap_bbox <- function(map, crs = 3857) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))

  # Coonvert the bbox to an sf polygon, transform it to 3857,
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(map_bbox, crs = 4326) %>%
    st_as_sfc() %>%
    st_transform(crs = crs) %>%
    st_bbox()

  # Overwrite the bbox of the ggmap object with the transformed coordinates
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}


library(ggmap)
get_map("Japan", zoom = 5, source = "stamen") %>% ggmap()

st_bbox(CCSBT %>% st_transform(crs = 4326))
# gm <- get_map(location = "New Zealand", zoom = 2, source = "google", maptype = "satellite")
gm <- get_map(location = c(lon = 0, lat = -30), zoom = 1, source = "google", maptype = "satellite")
ggmap(gm)


map <- ggmap_bbox(gm, crs = proj_ccsbt())
# map <- ggmap_bbox(gm)
# p <- ggmap(map) +
p <- ggplot() +
  annotation_map_tile(zoom = 7) +
  geom_sf(data = CCSBT, fill = NA)
  # coord_ccsbt()
p

@quantifish
Copy link
Owner Author

library(ncdf4)
library(ncdf.tools)
library(ecmwfr)

wf_set_key(user = "youremail@gmail.com", key = "secretkey", service = "webapi")

request <- list(dataset      = "reanalysis-era5-single-levels-monthly-means",
                product_type = "reanalysis",
                variable     = "sea_surface_temperature",
                year         = 1979:2020,
                month        = 1:12,
                area         = "-30/160/-55/185",
                format       = "netcdf",
                target       = "era5-nz_sst_to_2020.nc")

ncfile <- wf_request(user = "41999", request = request, transfer = TRUE, path = "", verbose = FALSE)



library(raster)
# should add code to get this direct from web
data(era5_nz_sst)
prj <- "+proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

suppressWarnings({
  sst1 <- era5_nz_sst %>%
    projectRaster(crs = prj) %>%
    rasterToPoints()
})

sst <- sst1 %>%
  data.frame() %>%
  pivot_longer(cols = c(-x, -y)) %>%
  separate(col = name, into = c("year", "month", "day", "hr", "min", "sec"), sep = "\\.") %>%
  mutate(year = gsub("X", "", year), value = value - 273.15) %>%
  filter(year %in% c("1979", "1989", "2009", "2019"), month %in% c("01"))
nz <- get_coast(proj = prj, resolution = "med")

ggplot() +
  geom_tile(data = sst, aes(x = x, y = y, fill = value), alpha = 0.95) +
  facet_wrap(year ~ month) +
  geom_sf(data = nz, fill = "black", colour = NA, size = 0.3) +
  coord_sf() +
  scale_fill_viridis("SST (°C)", alpha = 0.8, option = "plasma")

@quantifish
Copy link
Owner Author

f <- "../data-raw/mfe-average-seasurface-temperature-19932012-GTiff.zip"
fz <- unzip(zipfile = f, list = TRUE)
fz
unzip(zipfile = f)
mfe_average_sst <- raster::raster(x = "average-seasurface-temperature-19932012.tif", values = TRUE)# %>%
projectRaster(crs = proj_nzsf)
names(mfe_average_sst) <- "layer"
mfe_average_sst[mfe_average_sst[] < -10 | mfe_average_sst[] > 35] <- NA
range(mfe_average_sst[], na.rm = TRUE)
plot(mfe_average_sst)
file.remove(fz$Name)
use_data(mfe_average_sst, overwrite = TRUE)
st_transform(3994) %>%
data(mfe_average_sst)
rpts <- mfe_average_sst %>%
  rasterToPoints() %>%
  data.frame()
ggplot() +
  geom_raster(data = rpts, aes(x = x, y = y, fill = layer)) +
  # plot_statistical_areas(area = "EEZ", fill = NA) +
  # plot_coast(proj = projection(mfe_average_sst), resolution = "med", fill = "black", colour = NA, size = 0.3) +
  coord_sf() +
  scale_fill_viridis(alpha = 0.8, option = "magma") +
  labs(fill = "SST (°C)")

@quantifish
Copy link
Owner Author

data(gebco_depth_raster)

rpts <- gebco_depth_raster %>%
  rasterToPoints() %>%
  data.frame()

ggplot() +
  geom_raster(data = rpts, aes(x = x, y = y, fill = layer)) +
  plot_statistical_areas(area = "EEZ", fill = NA) +
  plot_coast(resolution = "med", fill = "black", colour = NA, size = 0.3) +
  coord_sf() +
  scale_fill_viridis(alpha = 0.8, option = "magma") +
  labs(fill = "Depth (m)")

@quantifish
Copy link
Owner Author

x <- gebco_CCAMLR %>%
  st_as_stars() %>%
  st_transform(crs = proj_ccamlr()) %>%
  as.data.frame()
as("raster") %>%
as("SpatialPixelsDataFrame")
x <- gebco_CCAMLR %>%
  as("SpatialPixelsDataFrame") %>%
  st_transform(crs = proj_ccamlr())

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant