Transform or convert coordinates of simple feature

Source: R/transform.R

Transform or convert coordinates of simple feature

st_transform(x, crs, ...)

# S3 method for sfc
st_transform(
  x,
  crs = st_crs(x),
  ...,
  aoi = numeric(0),
  pipeline = character(0),
  reverse = FALSE,
  partial = TRUE,
  check = FALSE
)

# S3 method for sf
st_transform(x, crs = st_crs(x), ...)

# S3 method for sfg
st_transform(x, crs = st_crs(x), ...)

sf_proj_info(type = "proj", path)

st_wrap_dateline(x, options, quiet)

# S3 method for sfc
st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE)

# S3 method for sf
st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE)

# S3 method for sfg
st_wrap_dateline(x, options = "WRAPDATELINE=YES", quiet = TRUE)

Arguments

x

object of class sf, sfc or sfg
crs

coordinate reference system: integer with the EPSG code, or character with proj4string
...

ignored
aoi

area of interest, in degrees: WestLongitude, SouthLatitude, EastLongitude, NorthLatitude
pipeline

character; proj4 or WKT coordinate operation, to override the default operation
reverse

boolean; if TRUE, the inverse operation of the pipeline is applied
partial

logical; allow for partial projection, if not all points of a geometry can be projected (corresponds to setting environment variable OGR_ENABLE_PARTIAL_REPROJECTION to TRUE)
check

logical; perform a sanity check on resulting polygons?
type

character; one of have_datum_files, proj, ellps, datum, units or prime_meridians; see Details.
path

character; PROJ search path to be set
options

character; should have "WRAPDATELINE=YES" to function; another parameter that is used is "DATELINEOFFSET=10" (where 10 is the default value)
quiet

logical; print options after they have been parsed?

Details

Transforms coordinates of object to new projection. Features that cannot be transformed are returned as empty geometries.

Projecting to projections not supported by GDAL may be done by st_transform_proj, part of package lwgeom.

The st_transform method for sfg objects assumes that the CRS of the object is available as an attribute of that name.

sf_proj_info lists the available projections, ellipses, datums, units, or data search path of the PROJ library when type is equal to proj, ellps, datum, units or path; when type equals have_datum_files a boolean is returned indicating whether datum files are installed and accessible (checking for conus).

for PROJ >= 6, sf_proj_info does not provide option type = "datums". PROJ < 6 does not provide the option type = "prime_meridians".

for PROJ >= 7.1.0, the "units" query of sf_proj_info returns the to_meter variable as numeric, previous versions return a character vector containing a numeric expression.

For a discussion of using options, see https://github.com/r-spatial/sf/issues/280 and https://github.com/r-spatial/sf/issues/541

Examples

p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
sfc

#> Geometry set for 2 features 
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -30 ymin: 20 xmax: 7 ymax: 52
#> CRS:            EPSG:4326
#> POINT (7 52)
#> POINT (-30 20)
st_transform(sfc, 3857)

#> Geometry set for 2 features 
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
#> CRS:            EPSG:3857
#> POINT (779236.4 6800125)
#> POINT (-3339585 2273031)
st_transform(st_sf(a=2:1, geom=sfc), "+init=epsg:3857")

#> Simple feature collection with 2 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
#> CRS:            +init=epsg:3857
#>   a                     geom
#> 1 2 POINT (779236.4 6800125)
#> 2 1 POINT (-3339585 2273031)
try(st_transform(sfc, 3857, aoi = c(-280,-90,180,90)))

#> Error in CPL_transform(x, crs, aoi, pipeline, reverse) : 
#>   pipeline or area of interest require GDAL >= 3
if (sf_extSoftVersion()["GDAL"] >= "3.0.0") {
  st_transform(sfc, pipeline =
    "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes
  st_transform(sfc, pipeline =
    "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes
}
nc = st_read(system.file("shape/nc.shp", package="sf"))

#> Reading layer `nc' from data source `/tmp/RtmpKaoQvj/temp_libpath71c86f876e69/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> CRS:            4267
st_area(nc[1,]) # area from long/lat

#> 1137388604 [m^2]
st_area(st_transform(nc[1,], 32119)) # NC state plane, m

#> 1137598162 [m^2]
st_area(st_transform(nc[1,], 2264)) # NC state plane, US foot

#> 12244955726 [US_survey_foot^2]
library(units)

#> udunits system database from /usr/share/xml/udunits
set_units(st_area(st_transform(nc[1,], 2264)), m^2)

#> 1137598162 [m^2]
st_transform(structure(p1, proj4string = "+init=epsg:4326"), "+init=epsg:3857")

#> POINT (779236.4 6800125)
sf_proj_info("datum")

#>             name   ellipse
#> 1          WGS84     WGS84
#> 2         GGRS87     GRS80
#> 3          NAD83     GRS80
#> 4          NAD27    clrk66
#> 5        potsdam    bessel
#> 6       carthage clrk80ign
#> 7  hermannskogel    bessel
#> 8          ire65  mod_airy
#> 9         nzgd49      intl
#> 10        OSGB36      airy
#>                                                        definition
#> 1                                                   towgs84=0,0,0
#> 2                                    towgs84=-199.87,74.79,246.62
#> 3                                                   towgs84=0,0,0
#> 4               nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
#> 5                 towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
#> 6                                        towgs84=-263.0,6.0,431.0
#> 7         towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232
#> 8      towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
#> 9              towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
#> 10 towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894
#>                             description
#> 1                                      
#> 2  Greek_Geodetic_Reference_System_1987
#> 3             North_American_Datum_1983
#> 4             North_American_Datum_1927
#> 5           Potsdam Rauenberg 1950 DHDN
#> 6                 Carthage 1934 Tunisia
#> 7                         Hermannskogel
#> 8                          Ireland 1965
#> 9       New Zealand Geodetic Datum 1949
#> 10                            Airy 1830
st_wrap_dateline(st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326))

#> Geometry set for 1 feature 
#> geometry type:  MULTILINESTRING
#> dimension:      XY
#> bbox:           xmin: -180 ymin: 0 xmax: 180 ymax: 0
#> CRS:            EPSG:4326
#> MULTILINESTRING ((-179 0, -180 0), (180 0, 179 0))
library(maps)
wrld <- st_as_sf(maps::map("world", fill = TRUE, plot = FALSE))
wrld_wrap <- st_wrap_dateline(wrld, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"),
   quiet = TRUE)

#> Warning: GDAL Error 1: IllegalArgumentException: Points of LinearRing do not form a closed linestring
#> Warning: GDAL Error 1: IllegalArgumentException: Points of LinearRing do not form a closed linestring
wrld_moll <- st_transform(wrld_wrap, "+proj=moll")
plot(st_geometry(wrld_moll), col = "transparent")