From 6f600332c377f783533b1be5a4c355fb3dd9031a Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Tue, 16 Jul 2024 05:08:19 -0700 Subject: [PATCH] Dev (#95) * small fixes for Region = [shapefile] * Update convert_shapefile.R * Update DESCRIPTION --------- Co-authored-by: Jim Thorson --- DESCRIPTION | 4 +-- R/convert_shapefile.R | 53 +++++++++++++++++++++---------------- R/make_extrapolation_info.R | 14 +++++----- 3 files changed, 38 insertions(+), 33 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index aa41a84..ee3f21e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: FishStatsUtils Type: Package Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox -Version: 2.13.0 -Date: 2024-01-10 +Version: 2.13.1 +Date: 2024-07-15 Authors@R: c(person(given = "James", family = "Thorson", diff --git a/R/convert_shapefile.R b/R/convert_shapefile.R index 0a1f1ac..f7cdaef 100644 --- a/R/convert_shapefile.R +++ b/R/convert_shapefile.R @@ -32,51 +32,58 @@ convert_shapefile = function( file_path, #shapefile_input = rgdal::readOGR( file_path, verbose=FALSE, p4s=projargs_for_shapefile ) shapefile_input = sf::st_read( file_path ) - st_crs(shapefile_input) = st_crs(projargs_for_shapefile) + + # Need to run sf::st_make_valid to ensure sf::st_area works later + shapefile_input = sf::st_make_valid( shapefile_input ) + # sf::st_crs(shapefile_input) = sf::st_crs(projargs_for_shapefile) # rgdal::writeOGR - message("Reading shapefile with projargs: ", print(shapefile_input@proj4string)) - # raster::shapefile(.) has simplified read-write interface for future reference - if( !(shapefile_input@class %in% c("SpatialPolygonsDataFrame","SpatialPolygons")) ){ - warning( "object at `file_path` doesn't appear to be a shapefile") + message("Reading shapefile with projargs: ", print(sf::st_crs(shapefile_input))) + # + if( is.null(sf::st_crs(shapefile_input)) ){ + sf::st_crs(shapefile_input) = sf::st_crs(projargs_for_shapefile) } + # raster::shapefile(.) has simplified read-write interface for future reference + #if( !(shapefile_input %in% c("SpatialPolygonsDataFrame","SpatialPolygons")) ){ + # warning( "object at `file_path` doesn't appear to be a shapefile") + #} proj_orig = "+proj=longlat +ellps=WGS84 +no_defs" - shapefile_orig = sp::spTransform(shapefile_input, CRSobj=sp::CRS(proj_orig) ) + shapefile_orig = sf::st_transform(shapefile_input, crs=sf::st_crs(proj_orig) ) #shapefile_orig@proj4string = sp::CRS(proj_orig) # Infer projargs if missing, and project - utm_zone = floor((mean(sp::bbox(shapefile_orig)[1,]) + 180) / 6) + 1 + utm_zone = floor((mean(sf::st_bbox(shapefile_orig)[c('xmin','xmax')]) + 180) / 6) + 1 projargs_utm = paste0("+proj=utm +zone=",utm_zone," +ellps=WGS84 +datum=WGS84 +units=km +no_defs ") if( is.null(projargs) || is.na(projargs) ){ projargs = projargs_utm } - shapefile_proj = sp::spTransform(shapefile_orig, CRSobj=sp::CRS(projargs) ) + shapefile_proj = sf::st_transform(shapefile_orig, crs=sf::st_crs(projargs) ) # Check for total area - Total_area = sum(raster::area(shapefile_proj)) - if( Total_area < 100 ){ + Total_area = as.numeric(sum(sf::st_area(shapefile_proj))) + if( (Total_area) < 100 ){ warning("Total area for polygons in shapefile is ", Total_area, " km^2 and this might indicate an issue with identifying or specifying the projection" ) } # Determine bounds for box - bbox = shapefile_proj@bbox - bbox[,1] = floor(bbox[,1] / grid_dim_km[1]) * grid_dim_km[1] - bbox[,2] = ceiling(bbox[,2] / grid_dim_km[2]) * grid_dim_km[2] + bbox = sf::st_bbox(shapefile_proj) + bbox[c('xmin','xmax')] = floor(bbox[c('xmin','xmax')] / grid_dim_km[1]) * grid_dim_km[1] + bbox[c('ymin','ymax')] = ceiling(bbox[c('ymin','ymax')] / grid_dim_km[2]) * grid_dim_km[2] # Make grid - grid_proj = expand.grid(X=seq(bbox[1,1],bbox[1,2],by=grid_dim_km[1]), Y=seq(bbox[2,1],bbox[2,2],by=grid_dim_km[2])) - grid_proj = sp::SpatialPointsDataFrame( coords=grid_proj[,c("X","Y")], data=grid_proj[,1:2], proj4string=sp::CRS(projargs) ) + grid_proj = expand.grid(X=seq(bbox['xmin'],bbox['xmax'],by=grid_dim_km[1]), Y=seq(bbox['ymin'],bbox['ymax'],by=grid_dim_km[2])) + grid_proj = sf::st_as_sf( grid_proj[,c("X","Y")], coords=1:2, crs=sf::st_crs(projargs) ) # Restrict points to those within a polygon - grid_proj@data = sp::over(grid_proj, shapefile_proj) - grid_proj = subset( grid_proj, !apply(grid_proj@data,MARGIN=1,FUN=function(vec){all(is.na(vec))}) ) + grid_proj = sf::st_intersection(grid_proj, shapefile_proj) + grid_proj_coordinates = sf::st_coordinates(grid_proj) # Return to original coordinates - grid_orig = sp::spTransform( grid_proj, CRSobj=sp::CRS(proj_orig) ) - grid_orig = as.data.frame(grid_orig) + grid_orig = sf::st_transform( grid_proj, crs=sf::st_crs(proj_orig) ) + grid_orig = sf::st_coordinates(grid_orig) # Combine extrapolation_grid = data.frame( "Lat"=grid_orig[,'Y'], "Lon"=grid_orig[,'X'], "Area_km2"=prod(grid_dim_km), # "Stratum"=grid_orig[,'AreaName'], - "Include"=1, "E_km"=grid_proj@coords[,'X'], "N_km"=grid_proj@coords[,'Y'] ) + "Include"=1, "E_km"=grid_proj_coordinates[,'X'], "N_km"=grid_proj_coordinates[,'Y'] ) # make plots if( make_plots==TRUE ){ @@ -100,9 +107,9 @@ convert_shapefile = function( file_path, } # Compare areas - area_shapefile_input = sum( raster::area(shapefile_input) ) - area_shapefile_orig = sum( raster::area(shapefile_orig) ) / 1000^2 # Convert to square-kiometers - area_shapefile_proj = sum( raster::area(shapefile_proj) ) + area_shapefile_input = as.numeric(sum( sf::st_area(shapefile_input) )) + area_shapefile_orig = as.numeric(sum( sf::st_area(shapefile_orig) )) / 1000^2 # Convert to square-kiometers + area_shapefile_proj = as.numeric(sum( sf::st_area(shapefile_proj) )) area_grid_proj = sum( extrapolation_grid[,'Area_km2'] ) # Messages diff --git a/R/make_extrapolation_info.R b/R/make_extrapolation_info.R index 3700a57..41ac2ff 100644 --- a/R/make_extrapolation_info.R +++ b/R/make_extrapolation_info.R @@ -305,21 +305,19 @@ make_extrapolation_info = function( Region, plot.make_extrapolation_info <- function(x, cex=0.01, land_color="grey", map_resolution="medium", ...) { par( mfrow=c(1,2), mar=c(3,3,2,0), mgp=c(1.75,0.25,0) ) + require(sf) # CRS for original and new projections - CRS_orig = sp::CRS( '+proj=longlat' ) - CRS_proj = sp::CRS( x$projargs ) + CRS_orig = st_crs( '+proj=longlat' ) + CRS_proj = st_crs( x$projargs ) # Data for mapping - map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50 )) - # Fix warning messages from projecting rnaturalearth object - # Solution: Recreate SpatialPolygonsDataFrame from output - map_data = sp::SpatialPolygonsDataFrame( Sr=sp::SpatialPolygons(slot(map_data,"polygons"),proj4string=CRS_orig), data=slot(map_data,"data") ) + map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50 ), return="sf") # Plot #1 -- Latitude plot( x$Data_Extrap[which(strip_units(x$Area_km2_x)>0),c('Lon','Lat')], cex=cex, main="Extrapolation (Lat-Lon)", ... ) - map_data_orig = sp::spTransform(map_data, CRSobj=CRS_orig) - sp::plot( map_data_orig, col=land_color, add=TRUE ) + map_data_orig = st_transform(map_data, crs=CRS_orig) + plot( map_data_orig, col=land_color, add=TRUE ) # Plot #2 -- Projection coordinates if( !any(is.na(x$Data_Extrap[,c('E_km','N_km')])) ){