DOI-USGS / intersectr

See official repository here: https://code.usgs.gov/water/intersectr
Creative Commons Zero v1.0 Universal
7 stars 5 forks source link

Example disparate points against metdata. #13

Closed dblodgett-usgs closed 5 years ago

dblodgett-usgs commented 5 years ago
library(sf)
library(dplyr)
library(intersectr)
library(ncmeta)
library(RNetCDF)
p <- st_sfc(st_point(c(-93.2, 47.2)), 
            st_point(c(-86.7, 42.4)), 
            crs = 4326)

geom <- st_sf(id = c("point1", "point2"), geom = p, stringsAsFactors = FALSE) %>%
  st_transform(5070) %>%
  st_buffer(dist = 20)

nc_file <- "https://cida.usgs.gov/thredds/dodsC/UofIMETDATA"
(nc_var <- nc_vars(nc_file))

variable_name <- "precipitation_amount"
(nc_coord_vars <- nc_coord_var(nc_file, variable_name))

(nc_prj <- nc_gm_to_prj(nc_grid_mapping_atts(nc_file)))

nc <- open.nc(nc_file)
X_coords <- var.get.nc(nc, nc_coord_vars$X, unpack = TRUE)
Y_coords <- var.get.nc(nc, nc_coord_vars$Y, unpack = TRUE)

start_date <- "2010-01-01 00:00:00"
end_date <- "2010-01-10 00:00:00"

intersected <- list()

for(point in 1:nrow(geom)) {
  geom_step <- geom[point, ]
  (cell_geometry <-
     create_cell_geometry(X_coords = X_coords,
                          Y_coords = Y_coords,
                          prj = nc_prj,
                          geom = st_transform(geom_step, nc_prj),
                          buffer_dist = 1000,
                          regularize = TRUE))

  data_source_cells <- st_sf(select(cell_geometry, grid_ids))
  target_polygons <- st_sf(select(geom_step, id))
  st_agr(data_source_cells) <- "constant"
  st_agr(target_polygons) <- "constant"
  area_weights <- calculate_area_intersection_weights(
    data_source_cells,
    target_polygons)

  intersected[[point]] <- execute_intersection(nc_file = nc_file,
                                               variable_name = variable_name,
                                               intersection_weights = area_weights,
                                               cell_geometry = cell_geometry, 
                                               x_var = nc_coord_vars$X,
                                               y_var = nc_coord_vars$Y,
                                               t_var = nc_coord_vars$T, 
                                               start_datetime = start_date, 
                                               end_datetime = end_date, 
                                               status = TRUE)

}