## read all shapefiles, we
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
library(tabularaster)
## tidync@expt-group-by
devtools::load_all(".")
## Loading tidync
bom <- "/rdsi/public/CFA/BoM_daily_vars/tmax/bom-tmax_day-19610101-19610131.nc"
tidync(bom) %>% activate(tmax_day)
## Variables: tmax_day, (crs, lat_bounds, lon_bounds, time_bounds)
## Dimensions:
## # A tibble: 3 x 5
## variable_name .variable_ .dimension_ dimension_name dimension_length
## <chr> <dbl> <int> <chr> <int>
## 1 tmax_day 7 2 longitude 841
## 2 tmax_day 7 0 latitude 681
## 3 tmax_day 7 3 time 31
shp <-read_sf("/mnt/WineAustralia_working_files/shapes/Wine-Regions-GI-Regions/individual_shp/region_tasmania.shp")
library(ggplot2)
## this point-in-polygon test needs review here, sp is uncharacteristically slow
system.time(d <- tidync(bom) %>% activate(tmax_day) %>% hyper_group_by(shp) %>% hyper_tibble())
## Warning: projections not the same
## x: NA
## query: +proj=longlat +ellps=GRS80 +no_defs
## user system elapsed
## 14.444 1.380 15.860
d
## # A tibble: 91,636 x 4
## tmax_day longitude latitude time
## <dbl> <dbl> <dbl> <dbl>
## 1 21.20 147.30 -39.45 58804
## 2 18.25 143.95 -39.60 58804
## 3 18.24 143.95 -39.65 58804
## 4 18.32 144.00 -39.65 58804
## 5 18.30 144.05 -39.65 58804
## 6 21.02 148.00 -39.65 58804
## 7 18.31 143.90 -39.70 58804
## 8 18.25 143.95 -39.70 58804
## 9 18.30 144.00 -39.70 58804
## 10 18.34 144.05 -39.70 58804
## # ... with 91,626 more rows
ggplot(d, aes(longitude, latitude, fill = tmax_day)) + geom_raster() + facet_wrap(~time) + coord_equal()
shp <-read_sf("/mnt/WineAustralia_working_files/shapes/Wine-Regions-GI-Regions/individual_shp/region_adelaide_hills.shp")
system.time(d <- tidync(bom) %>% activate(tmax_day) %>% hyper_group_by(shp) %>% hyper_tibble())
## Warning: projections not the same
## x: NA
## query: +proj=longlat +ellps=GRS80 +no_defs
## user system elapsed
## 1.020 0.024 1.049
d
## # A tibble: 1,736 x 4
## tmax_day longitude latitude time
## <dbl> <dbl> <dbl> <dbl>
## 1 24.40 138.85 -34.70 58804
## 2 23.50 138.95 -34.70 58804
## 3 23.70 138.80 -34.75 58804
## 4 23.98 138.85 -34.75 58804
## 5 23.44 138.90 -34.75 58804
## 6 23.20 138.95 -34.75 58804
## 7 23.12 139.00 -34.75 58804
## 8 22.95 139.05 -34.75 58804
## 9 23.49 138.80 -34.80 58804
## 10 23.51 138.85 -34.80 58804
## # ... with 1,726 more rows
ggplot(d, aes(longitude, latitude, fill = tmax_day)) + geom_raster() + facet_wrap(~time)
