Reproject a vector "R" object or file

Reproject a vector "R" object or file to a different reference system. The function is a simple wrapper around sf::st_tranform with additional checks on inputs allowing to specify the output projection in several ways: 1: passing a valid proj4 string (e.g., reproj_vect(in_vect, "+init=epsg:4326") 2: passing a numeric or character that can be interpreted as an EPSG code (e.g., reproj_vect(in_vect, 4326)); 3: passing the name of a valid R vector or raster object: (e.g.,reproj_vect(in_vect, rast_obj), with rast_obj being an existing R object; 4: passing the path to a valid vector or raster file: EPSG code (e.g.,reproj_vect(in_vect, "D:/Temp/myfile.tif")

The reprojected vector can be automatically saved to a shapefile through write_shape by specifying a valid path with the out_file argument.

reproj_vect(in_vect, in_projobj, out_file = NULL,
  out_type = "vectobject", out_class = NULL, overwrite = FALSE,
  verbose = TRUE)

Arguments

in_vect	A vector object (*sf or sp), or the path to a vector file
in_projobj	R object or filename from which the output projection should be derived (see @description )
out_file	character Path where the reprojected vector should be saved. If NULL, the reprojected vector is not saved, but just sent back to the caller, Default: NULL
out_type	character ["vectfile" | "vectobject"] If "vectfile", and out_file is not NULL, the function returns the name of the saved shapefile. Otherwise, it returns the reprojected vector object, with the format specified by out_class, Default: 'vectobject'
out_class	character ["sf" | "sp"] Specifies if the reprojected object should have class sf or sp. If NULL, the returned reprojected object has the same class of in_vect, Default: NULL
overwrite	logical If TRUE, overwrite existing files, Default: FALSE
verbose	logical If FALSE, suppress processing messages, Default: TRUE

Value

a vector of class *sf or sp (depending on out_class), or the path of the file where the reprojected input was saved (if out_type == "vectfile")

See also

sf::st_write sf::st_transform write_shape

Examples

library(sprawl.data)
library(gridExtra)
#> 
#> Attaching package: ‘gridExtra’
#> The following object is masked from ‘package:dplyr’:
#> 
#>     combine
# reproject a vector file
in_vect <- system.file("extdata/shapes","lc_polys.shp",
                                   package = "sprawl.data")
in_vect <- read_vect(in_vect)
bounds  <- get_boundaries("PHL", 1)
#> get_boundaries --> Downloading data for: "PHL", Level: 1

message("Input projection is: ", get_proj4string(in_vect))
#> Input projection is: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
plot_vect(in_vect, fill_var = "category", borders_layer = bounds)
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.

# reproject to 3857 (web mercator)
out_proj <-  3857
out_vect <- reproj_vect(in_vect, 3857)
#> reproj_vect --> Reprojecting in_vect to +init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

message("Output projection is: ", get_proj4string(out_vect))
#> Output projection is: +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +nadgrids=@null +units=m +no_defs

# Do the same, but also save the output to file
out_file <- tempfile(fileext = ".shp")
out_proj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
out_vect <- reproj_vect(in_vect, 3857, out_file = out_file)
#> reproj_vect --> Reprojecting in_vect to +init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
read_vect(out_file)
#> Simple feature collection with 13 features and 4 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 13378130 ymin: 1706901 xmax: 13507440 ymax: 1815546
#> epsg (SRID):    NA
#> proj4string:    +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs
#> # A tibble: 13 x 5
#>       id lc_type   category sup_catego                                  geometry
#>    <dbl> <chr>     <chr>    <chr>                                  <POLYGON [m]>
#>  1     1 forest_1  Forest   Vegetation ((13439627 1779722, 13440374 1779722, 13…
#>  2     2 forest_2  Forest   Vegetation ((13399126 1735311, 13399217 1735673, 13…
#>  3     3 urban_1   Other    Other      ((13424474 1746623, 13424535 1746707, 13…
#>  4     4 cropland… Cropland Vegetation ((13429260 1742926, 13429620 1744112, 13…
#>  5     5 cropland… Cropland Vegetation ((13438831 1751932, 13440785 1752637, 13…
#>  6     6 forest_3  Forest   Vegetation ((13494772 1805481, 13496236 1806264, 13…
#>  7     7 reiver_b… Other    Other      ((13405185 1727398, 13405679 1727635, 13…
#>  8     8 cropland… Cropland Vegetation ((13444080 1731748, 13444335 1732055, 13…
#>  9     9 urban_2_… Other    Other      ((13423271 1708906, 13425121 1709004, 13…
#> 10    10 sparsfor… Forest   Vegetation ((13380132 1709289, 13381309 1709303, 13…
#> 11    11 water     Other    Other      ((13484992 1788081, 13487983 1786824, 13…
#> 12    12 small_fe… Other    Other      ((13445172 1797572, 13445287 1797560, 13…
#> 13    13 out_feat  Other    Other      ((13433826 1815546, 13435617 1814804, 13…

# use a different spatial file or object to set the output projection:
rast <- read_rast(system.file("extdata/MODIS_test", "EVIts_test.tif",
  package = "sprawl.data"))
out_vect <- reproj_vect(in_vect, rast)
#> reproj_vect --> Reprojecting in_vect to projection ofrast

message("Output projection is: ", get_proj4string(out_vect))
#> Output projection is: +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs

plot_vect(in_vect, borders_layer = bounds, fill_var = "category",
   title = "Original (lat/lon)")
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.

plot_vect(out_vect, borders_layer = bounds, fill_var = "category",
   title = "Reprojected (sinusoidal)")
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.