library(angstroms)
## Loading required package: raster
## Loading required package: sp
u <- "http://dapds00.nci.org.au/thredds/dodsC/rr6/eReefs/4kmReanalyses/ocean_avg_20060102.nc"
(vars <- ncmeta::nc_vars(u))
## # A tibble: 23 x 5
## id name type ndims natts
## <int> <chr> <chr> <int> <int>
## 1 0 s_rho NC_DOUBLE 1 7
## 2 1 ocean_time NC_DOUBLE 1 4
## 3 2 zeta NC_FLOAT 3 6
## 4 3 zeta_detided NC_FLOAT 3 6
## 5 4 ubar NC_FLOAT 3 6
## 6 5 vbar NC_FLOAT 3 6
## 7 6 u NC_FLOAT 4 7
## 8 7 v NC_FLOAT 4 7
## 9 8 temp NC_FLOAT 4 6
## 10 9 salt NC_FLOAT 4 5
## # … with 13 more rows
vars$name
## [1] "s_rho" "ocean_time" "zeta" "zeta_detided" "ubar"
## [6] "vbar" "u" "v" "temp" "salt"
## [11] "Hsbl" "Uwind" "Vwind" "shflux" "lat_rho"
## [16] "lat_u" "lat_v" "lon_rho" "lon_u" "lon_v"
## [21] "Cs_r" "angle" "h"
coords <- romscoords(u, c("lon_u", "lat_u"))
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named xi_u BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named eta_u BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named xi_u BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named eta_u BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## get a bit mappy
plot(values(coords), pch = ".")
maps::map(add = TRUE) ## look ok?
temp <- roms_xy(u, varname = "temp", slice = c(1, 1))
## ok let's plot properly (things have improved, most of these ideas came from angstroms ...)
library(anglr)
## for this to work we have to pretend that lon_u/lat_u are the mass points, they are close enough if we drop one column
mesh_plot(crop(temp, extent(temp, 1, nrow(temp), 1, ncol(temp)-1)), coords = coords)
maps::map(add = TRUE)
library(sp)
## now, let's draw a polygon in xy
#p <- raster::drawPoly() ## you might read a shapefile, whatever
m <- structure(c(151.787077738327, 153.276303417153, 154.590326074941,
150.735859612097, 149.290434688531, 151.787077738327, -18.3248255134419,
-18.5942335261058, -21.4037742296007, -21.0189056400808, -19.0560758335296,
-18.3248255134419), .Dim = c(6L, 2L))
p <- sp::SpatialPolygonsDataFrame(sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(m)), "1")), proj4string = sp::CRS("+proj=longlat")), data.frame(x = 1))
p0 <- romsmap(p, coords)
plot(p0) ## see how this is in index space now
plot(temp, add = TRUE, col = hcl.colors(10, alpha = 0.3))
## so, we can extract in index space
str(extract(temp, p0))
## List of 1
## $ : num [1:5759] 14.3 13.5 15 14.1 13.4 ...
## tada, etc. - there are analogous functions romsdata_xt, etc and they work IIRC but I haven't done this in a while
n <- 180000
xy <- cbind(runif(n, 145, 160), runif(n, -28, -8))
## it has to be an sp object sadly (but wouldn't be hard to fix this)
temp0 <- extract(temp, romsmap(sp::SpatialPointsDataFrame(sp::SpatialPoints(xy), data.frame(n = 1:n)), coords),
method = "bilinear")
## so for extract we have index, and we can mesh plot or plot in index space
## for our proper mapping we have geography, so we can flip between these as we like
plot(xy[!is.na(temp0), ], pch = 19, cex = 0.1, col = colourvalues::colour_values(temp0[!is.na(temp0)]),
asp = 1/cos(-18 * pi/180))
plot(p, add = TRUE, border = "hotpink", lwd = 4)
maps::map(add = TRUE)
## HTH