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

Modularizing a bit and adding a function to generate OPeNDAP urls. #16

Closed dblodgett-usgs closed 5 years ago

dblodgett-usgs commented 5 years ago

We can now do:

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)

(cell_geometry <-
    create_cell_geometry(X_coords = X_coords,
                         Y_coords = Y_coords,
                         prj = nc_prj,
                         geom = st_transform(geom, nc_prj),
                         regularize = TRUE))

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

dap <- intersectr:::get_dap_url(min(cell_geometry$X_ind), max(cell_geometry$X_ind),
                                min(cell_geometry$Y_ind), max(cell_geometry$Y_ind),
                                nc_file, variable_name, x_var, y_var, t_var,
                                start_date, end_date)

And...

$ nccopy https://cida.usgs.gov/thredds/dodsC/UofIMETDATA?lon[759:914],lat[417:532],day[11324:11326],precipitation_amount[11324:11326][417:532][759:914] test.nc