Execute intersection failing with CFSR data and Alaska HUC12s

[dblodgett-usgs]

Getting:

Error in seq.default(min(cell_geometry$X_ind), max(cell_geometry$X_ind),  : 
  'from' must be a finite number 
4. stop("'from' must be a finite number") 
3. seq.default(min(cell_geometry$X_ind), max(cell_geometry$X_ind), 1) 
2. seq(min(cell_geometry$X_ind), max(cell_geometry$X_ind), 1) at execute_intersection.R#117
1. execute_intersection(f, "A_PCP_L1_Accum_1", weights, cells, "lon",  "lat", "time") 

library(RNetCDF)
library(sf)
library(dplyr)
library(intersectr)

download.file("https://github.com/USGS-R/intersectr/files/3118725/CFSR_HUC12_AL.zip", 
              "CFSR_HUC12_AL.zip")
unzip("CFSR_HUC12_AL.zip")
f <- "pgbh06.gdas.19790101-19790105.grb2.nc"

ncmeta::nc_coord_var(f)

nc <- open.nc(f)
x <- var.get.nc(nc, "lon")
y <- var.get.nc(nc, "lat")
close.nc(nc)

prj <-  ncmeta::nc_gm_to_prj(ncmeta::nc_grid_mapping_atts(f))

geom <- read_sf("geom.gpkg") %>% 
  st_transform(prj)

cells <- create_cell_geometry(x, y,  prj, geom, regularize = TRUE)

weights <- calculate_area_intersection_weights(select(cells, grid_ids),
                                               select(geom, HUC12))

execute_intersection(f, "A_PCP_L1_Accum_1", weights, cells, "lon", "lat", "time") 

CFSR_HUC12_AL.zip

[dblodgett-usgs]

Traced the issue back to the create_cell_geometry step. X_inds <- which(X_coords > req_bbox$xmin & X_coords < req_bbox$xmax) is giving all kinds of false positives for lat/lon data that crosses the international date line.

[dblodgett-usgs]

Fixed in #12

[dblodgett-usgs]

Final script used to debug with some added visualization goodness.

library(RNetCDF)
library(sf)
library(dplyr)
library(intersectr)

download.file("https://github.com/USGS-R/intersectr/files/3118725/CFSR_HUC12_AL.zip", 
              "CFSR_HUC12_AL.zip")
unzip("CFSR_HUC12_AL.zip")
f <- "pgbh06.gdas.19790101-19790105.grb2.nc"

ncmeta::nc_coord_var(f)

nc <- open.nc(f)
x <- var.get.nc(nc, "lon")
y <- var.get.nc(nc, "lat")
close.nc(nc)

prj <-  ncmeta::nc_gm_to_prj(ncmeta::nc_grid_mapping_atts(f))

geom <- read_sf("SW_alaska_huc4/SW_alaska_huc4.shp") %>% 
  st_transform(prj)

cells <- create_cell_geometry(x, y,  prj, geom, buffer_dist = 0.5, regularize = TRUE)

cellsproj <- st_transform(cells, "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
geom <- st_transform(geom, st_crs(cellsproj)) %>%
  lwgeom::st_make_valid() %>%
  st_buffer(0)

p_geom <- st_simplify(geom, dTolerance = 50)

plot(cellsproj$geometry)
plot(p_geom$geometry, add = TRUE)

weights <- calculate_area_intersection_weights(select(cellsproj, grid_ids),
                                               select(geom, HUC12))

intersected <- execute_intersection(f, "A_PCP_L1_Accum_1", weights, cells, "lon", "lat", "time") 

geomHUCs <- geom$HUC12
intersectedHUCs <- colnames(intersected[,2:ncol(intersected)])

nonintersectedHUCs <- geomHUCs[which(!geomHUCs %in% intersectedHUCs)]

geomnotintersected <- geom[geom$HUC12 %in% nonintersectedHUCs,]

plot(st_geometry(geom), col = 'black', axes = TRUE)
plot(st_geometry(geomnotintersected), pch = 3, col = 'red', add = TRUE)