James-Thorson-NOAA / FishStatsUtils

Shared resources for spatio-temporal models
GNU General Public License v3.0
10 stars 16 forks source link

combine_extrapolation_info doesn't create strata_per_region with user-defined grid #71

Open Jason-Conner-NOAA opened 3 years ago

Jason-Conner-NOAA commented 3 years ago

Attempting to use combine_extrapolation_info() to combine 2 user defined grids for the EBS/NBS. `# EBS Grid prep EBSgrid <- read.csv(here::here("data","EBSThorsonGrid.csv")) EBS_input <- cbind(Lat=EBSgrid$Lat,Lon=EBSgrid$Lon,Area_km2=EBSgrid$Shape_Area/1000000) EBS_Extrapolation_List = make_extrapolation_info( Region=Region, input_grid=EBS_input, strata.limits=strata.limits, create_strata_per_region=FALSE, max_cells = max_cells )

# NBS grid prep
NBSgrid <- read.csv(here::here("data","NBSThorsonGrid.csv"))
NBS_input <- cbind(Lat=NBSgrid$Lat,Lon=NBSgrid$Lon,Area_km2=NBSgrid$Shape_Area/1000000)
NBS_Extrapolation_List = make_extrapolation_info( Region=Region, input_grid=NBS_input, strata.limits=strata.limits,
                                                  create_strata_per_region=FALSE, max_cells = max_cells )

# Combine grids
Extrapolation_List <- combine_extrapolation_info(c(EBS_Extrapolation_List, NBS_Extrapolation_List), create_strata_per_region=TRUE)`

This results in combining the areas, but does not separate the region subtotals in the Table_for_SS3 output.

James-Thorson-NOAA commented 3 years ago

Thanks for asking here where our discussion is searchable!

Please note first that the intended syntax is:

# Combine grids
Extrapolation_List <- combine_extrapolation_info( EBS_Extrapolation_List, NBS_Extrapolation_List, create_strata_per_region=TRUE)`

i.e., passing both extrapolation-grid objects via the ... argument. Also, please note that the user will then have to do additional work to make a combined-stratum column of Extrapolation_List$a_gl. I'm closing for now, but feel free to keep commenting here if that's useful.

Jason-Conner-NOAA commented 3 years ago

Thanks @James-Thorson-NOAA , that worked for producing the subarea calculations in SS3!

However, now plot_results does not produce the predicted density map plots, returning this error:

 #Making plots of spatial predictions
 #plot_num 3: Plotting density maps (in log-space)
Error in sp::SpatialPointsDataFrame(coords = loc_g, data = data.frame(y = Y_gt[,  : 
  row.names of data and dimnames of coords do not match
In addition: Warning messages:
1: In spTransform(xSP, CRSobj, ...) :
  NULL source CRS comment, falling back to PROJ string
2: In spTransform(xSP, CRSobj, ...) :
  NULL source CRS comment, falling back to PROJ string
3: In par(Par) : "strata_names" is not a graphical parameter
4: In sp::SpatialPointsDataFrame(coords = loc_g, data = data.frame(y = Y_gt[,  :
  forming a SpatialPointsDataFrame based on maching IDs, not on record order. Use match.ID = FALSE to match on record order
James-Thorson-NOAA commented 3 years ago

what is outputted when you do lapply(Extrapolation_List, dim) and lapply(Extrapolation_List, dim)?

I'm trying to figure out if the problem arises from the code I provided above, or from efforts to add the combined-stratum index.

Jason-Conner-NOAA commented 3 years ago

lapply(Extrapolation_List, dim) is listed twice, did you mean a different function for one of those?

The output is:

lapply(Extrapolation_List, dim) $a_el [1] 4000 3

$Data_Extrap [1] 4000 5

$zone NULL

$flip_around_dateline NULL

$projargs NULL

$Area_km2_x NULL

James-Thorson-NOAA commented 3 years ago

sorry I meant lapply(Extrapolation_List, length)

Jason-Conner-NOAA commented 3 years ago

Gotcha, output for lapply(Extrapolation_List, length):

lapply(Extrapolation_List, length) $a_el [1] 12000

$Data_Extrap [1] 20000

$zone [1] 1

$flip_around_dateline [1] 1

$projargs [1] 1

$Area_km2_x [1] 4000

James-Thorson-NOAA commented 3 years ago

well, nothing obviously screwy there, so I guess its just some weird interaction affecting the plots (which I don't maintain with as high of scrutiny).

I guess maybe email me the custom RDS files you're combining, and I can try running those through my existing script for pollock EBS+NBS..?