Spatial data in R

Change map CRS

Mercator ESPG:3857 - most common USA projection in google maps

md.map+coord_sf(crs = st_crs(3857))

# problem - polar regions proportionally larger because to flatten stretches poles to make a flat, rectangular map
# Greenland has an area of about 836,300 square miles, compared to South America’s 6.9 million square miles

Selecting a projection - depends on spatial scope of map, goals
Do you want to preserve precise angles and directions between objects or preserve areas? Usually cannot do all.

md.map+coord_sf(crs = st_crs(2163))# ESPG 2163 - USA equal area 'Albers' - better to show states proportionally -optimized for USA

For smaller areas precision is more important than distortion

md.map+coord_sf(crs = st_crs(2248))# NAD83 Maryland - http://spatialreference.org/ref/epsg/nad83-maryland/

Assign CRS to shapefiles

md<-md.state%>%st_transform(2248)
md.cos<-md.cos%>%st_transform(2248)
st_crs(md) # units are in us ft

## Coordinate Reference System:
##   EPSG: 2248 
##   proj4string: "+proj=lcc +lat_1=39.45 +lat_2=38.3 +lat_0=37.66666666666666 +lon_0=-77 +x_0=399999.9998983998 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

Plot MD shapefile

md.map<-ggplot(md)+
  geom_sf()
md.map

Convert coordinates to sf object

If we have longitude and latitude for our sites, we can project these into a spatial format to work with

# lat and lon of seedbank sites
sites<-data.frame(site = c('Blackwater','Elmwood','Anderson'),
                  lat = c(38.437283, 38.16083, 38.23056),
                  lon = c(-76.215753, -75.7961, -75.774267))


# convert coordinates to a sf object
sites.sf <- st_as_sf(sites, coords = c("lon", "lat"), crs=4326)%>%
  st_transform(2248) # assign crs
glimpse(sites.sf) # converting to sf drops coordinates

## Observations: 3
## Variables: 2
## $ site     <fct> Blackwater, Elmwood, Anderson
## $ geometry <POINT [US_survey_foot]> POINT (1536969 281612.6), POINT (165…

Plot coordiantes on map

md.map+
  geom_sf(data=sites.sf, size=2.5, shape=21, fill='red')

#### Find coordiantes
To add labels, need to give coordinates of where to place

# find coordinates for each site
st_coordinates(sites.sf) 

##         X        Y
## 1 1536969 281612.6
## 2 1658494 182251.6
## 3 1664430 207727.2

# add coordinates as lat and lon to sites.sf 
sites.sf$lon<-st_coordinates(sites.sf)[,1] # add longitude
sites.sf$lat<-st_coordinates(sites.sf)[,2] # add latitude

glimpse(sites.sf) # examine sf object

## Observations: 3
## Variables: 4
## $ site     <fct> Blackwater, Elmwood, Anderson
## $ geometry <POINT [US_survey_foot]> POINT (1536969 281612.6), POINT (165…
## $ lon      <dbl> 1536969, 1658494, 1664430
## $ lat      <dbl> 281612.6, 182251.6, 207727.2

Add labels to map

longitude is x axis, latitude is y axis

md.map+
  geom_sf(data=sites.sf, size=2.5, shape=21, fill='red')+
  geom_label(data=sites.sf, aes(lon, lat, label=site))

library(ggrepel) # spread overlapping labels

md.map+
  geom_sf(data=sites.sf, size=2.5, shape=21, fill='red')

  geom_label_repel(data=sites.sf, aes(lon, lat, label=site))

## mapping: x = ~lon, y = ~lat, label = ~site 
## geom_label_repel: parse = FALSE, box.padding = 0.25, label.padding = 0.25, point.padding = 1e-06, label.r = 0.15, label.size = 0.25, segment.colour = NULL, segment.size = 0.5, segment.alpha = NULL, min.segment.length = 0.5, arrow = NULL, na.rm = FALSE, force = 1, max.iter = 2000, nudge_x = 0, nudge_y = 0, xlim = c(NA, NA), ylim = c(NA, NA), direction = both, seed = NA
## stat_identity: na.rm = FALSE
## position_identity

Calculating distances

st_distance() calculates the distance between points

st_distance(sites.sf) # all pairwise distances between sites

## Units: [US_survey_foot]
##          [,1]      [,2]      [,3]
## [1,]      0.0 156974.35 147327.79
## [2,] 156974.4      0.00  26158.06
## [3,] 147327.8  26158.06      0.00

Determine nearest distance to county boundary:

md.pt<-md%>%st_cast('POINT') # turn MD border into points
min(st_distance(md.pt, sites.sf)) # find the minimum distance and point along border 

## 5832.274 [US_survey_foot]

st_cast converts geometries to other types

In a similar manner we can add labels to each county polygon:

## find the centroid of a polygon 
# add county labels - calculate longitude and latitude for center of each county
glimpse(md.cos)

## Observations: 24
## Variables: 13
## $ statefp           <chr> "24", "24", "24", "24", "24", "24", "24", "24"…
## $ countyfp          <chr> "003", "009", "011", "023", "025", "027", "043…
## $ countyns          <chr> "01710958", "01676636", "00595737", "01711058"…
## $ affgeoid          <chr> "0500000US24003", "0500000US24009", "0500000US…
## $ geoid             <chr> "24003", "24009", "24011", "24023", "24025", "…
## $ name              <chr> "Anne Arundel", "Calvert", "Caroline", "Garret…
## $ lsad              <chr> "06", "06", "06", "06", "06", "06", "06", "06"…
## $ aland             <dbl> 1074553083, 552178304, 827350236, 1676036320, …
## $ awater            <dbl> 447833714, 341560888, 16777064, 22315761, 2319…
## $ state_name        <chr> "Maryland", "Maryland", "Maryland", "Maryland"…
## $ state_abbr        <chr> "MD", "MD", "MD", "MD", "MD", "MD", "MD", "MD"…
## $ jurisdiction_type <chr> "state", "state", "state", "state", "state", "…
## $ geometry          <MULTIPOLYGON [US_survey_foot]> MULTIPOLYGON (((1357…

# find centroid of each county
st_centroid(md.cos)

# use lon & lat for each centroid to place labels 
md.cos$lon<-st_coordinates(st_centroid(md.cos))[,1] # add longitude to sf
md.cos$lat<-st_coordinates(st_centroid(md.cos))[,2] # add latitude to sf

# create mape of MD counties colored by land area (aland)
md.co.map<-ggplot()+
  geom_sf(data=md, fill=NA, color='black')+
  geom_sf(data=md.cos, aes(fill=aland))+
  theme_bw()+
  scale_fill_viridis(alpha=.5)+
  theme(legend.position='none')
md.co.map

# label each county
md.co.map+geom_label_repel(data=md.cos, aes(x=lon, y=lat, label=name))

subsetting vector data

st_intersection finds the intersecting area between two vector objects

## determine which county each site is in, add to sf
sites.sf<-st_intersection(sites.sf, md.cos) # what counties do sites intersect with?
sites.sf$name # named counties

## [1] "Somerset"   "Somerset"   "Dorchester"

# filter shapefile to desired counties
site.cos<-md.cos%>%filter(name %in% sites.sf$name)

# crop map to just these counties
# recreate same map
co.map<-ggplot()+
  geom_sf(data=site.cos)+
  geom_sf(data=sites.sf, size=2, shape=21, fill='red', color='black')+
  geom_label_repel(data=sites.sf, aes(lon, lat, label=site))
co.map # counties missing wihtin map space - better to crop the map extent to show whole region

subsetting extent

st_crop clips vector data based on given limits (bounding box)

# determine full extent around sites
sites.bbox<-st_bbox(sites.sf)# bounding box
co.bbox<-st_crop(md.cos, sites.bbox)
glimpse(co.bbox) # US_survey_foot means spatial data is in US ft

## Observations: 3
## Variables: 15
## $ statefp           <chr> "24", "24", "24"
## $ countyfp          <chr> "039", "045", "019"
## $ countyns          <chr> "00596907", "01668606", "00596495"
## $ affgeoid          <chr> "0500000US24039", "0500000US24045", "0500000US…
## $ geoid             <chr> "24039", "24045", "24019"
## $ name              <chr> "Somerset", "Wicomico", "Dorchester"
## $ lsad              <chr> "06", "06", "06"
## $ aland             <dbl> 828090979, 969736590, 1400542033
## $ awater            <dbl> 752714478, 66797770, 1145384766
## $ state_name        <chr> "Maryland", "Maryland", "Maryland"
## $ state_abbr        <chr> "MD", "MD", "MD"
## $ jurisdiction_type <chr> "state", "state", "state"
## $ lon               <dbl> 1671442, 1707727, 1594945
## $ lat               <dbl> 165995.8, 260349.9, 298767.1
## $ geometry          <GEOMETRY [US_survey_foot]> POLYGON ((1613938 182251…

## zoom in map extent to sites
ggplot()+
  geom_sf(data=co.bbox)+
  geom_sf(data=sites.sf)

# points are on edge - need a buffer around to look nice
# repeat same steps but use st_buffer() to add space around bbox
sites.bbox<-st_bbox(st_buffer(sites.sf, 20000))# since map is in ft
co.bbox<-st_crop(md.cos, sites.bbox)

ggplot()+
  geom_sf(data=co.bbox, fill='darkgreen', alpha=.4)+
  geom_sf(data=sites.sf)+
  geom_text(data=sites.sf, aes(label=site, lon, lat), hjust=0.25, vjust=-1)+
  theme_bw()

Read in DEM raster

DEM = digital elevation model
data source: (http://data.imap.maryland.gov/datasets/1c6ce663eb3b499b9495010cb8c89df6v)

# download .tif file to working directory, read in
download.file('https://github.com/collnell/GWU-visual/blob/master/MD_DEM.tiff?raw=true', destfile=getwd())
dem<-raster('MD_DEM.tiff')

# look at file
dem

## class       : RasterLayer 
## dimensions  : 400, 400, 160000  (nrow, ncol, ncell)
## resolution  : 3252.954, 3252.954  (x, y)
## extent      : 585182.1, 1886364, -226561.9, 1074620  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=39.45 +lat_2=38.3 +lat_0=37.66666666666666 +lon_0=-77 +x_0=399999.9998983998 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs 
## data source : /Users/collnell/Dropbox/rstats/data viz/MD_DEM.tiff 
## names       : MD_DEM

dem@data

## An object of class ".SingleLayerData"
## Slot "values":
## logical(0)
## 
## Slot "offset":
## [1] 0
## 
## Slot "gain":
## [1] 1
## 
## Slot "inmemory":
## [1] FALSE
## 
## Slot "fromdisk":
## [1] TRUE
## 
## Slot "isfactor":
## [1] FALSE
## 
## Slot "attributes":
## list()
## 
## Slot "haveminmax":
## [1] FALSE
## 
## Slot "min":
## [1] Inf
## 
## Slot "max":
## [1] -Inf
## 
## Slot "band":
## [1] 1
## 
## Slot "unit":
## [1] ""
## 
## Slot "names":
## [1] "MD_DEM"

convert to dataframe to use with ggplot2

dem.df<-as.data.frame(dem, xy=TRUE)%>%filter(MD_DEM > 0)
str(dem.df)

## 'data.frame':    26858 obs. of  3 variables:
##  $ x     : num  648615 651868 655121 658374 661626 ...
##  $ y     : num  760710 760710 760710 760710 760710 ...
##  $ MD_DEM: num  580 538 621 600 695 ...

# plot raster values in histogram
ggplot(dem.df, aes(MD_DEM))+
  geom_histogram()

plot DEM raster pixel values

ggplot(dem.df)+
  geom_raster(data=dem.df, aes(x, y))

# assign color scale to raster values
ggplot(dem.df)+
  geom_raster(data=dem.df, aes(x, y, fill=MD_DEM), alpha=.8)+
  scale_fill_viridis(option='magma', begin=.4)

#### Crop raster to extent of sites

dem.bbox<-crop(dem,co.bbox)
dem.df<-as.data.frame(dem.bbox, xy=TRUE)%>%filter(MD_DEM != 0)

# plot it with sites
ggplot(dem.df)+
  geom_raster(data=dem.df, aes(x, y, fill=MD_DEM))+
  scale_fill_viridis(option='magma', begin=.4, 'DEM')+
  geom_sf(data=sites.sf, color='white', size=4)

extract DEM data for sites

extract pulls data values from raster layers that correspond to vector objects

sites.sf$DEM<-extract(dem.bbox, sites.sf, buffer=100, fun=mean)
glimpse(sites.sf)

## Observations: 3
## Variables: 19
## $ site              <fct> Elmwood, Anderson, Blackwater
## $ lon               <dbl> 1658494, 1664430, 1536969
## $ lat               <dbl> 182251.6, 207727.2, 281612.6
## $ statefp           <chr> "24", "24", "24"
## $ countyfp          <chr> "039", "039", "019"
## $ countyns          <chr> "00596907", "00596907", "00596495"
## $ affgeoid          <chr> "0500000US24039", "0500000US24039", "0500000US…
## $ geoid             <chr> "24039", "24039", "24019"
## $ name              <chr> "Somerset", "Somerset", "Dorchester"
## $ lsad              <chr> "06", "06", "06"
## $ aland             <dbl> 828090979, 828090979, 1400542033
## $ awater            <dbl> 752714478, 752714478, 1145384766
## $ state_name        <chr> "Maryland", "Maryland", "Maryland"
## $ state_abbr        <chr> "MD", "MD", "MD"
## $ jurisdiction_type <chr> "state", "state", "state"
## $ lon.1             <dbl> 1671442, 1671442, 1594945
## $ lat.1             <dbl> 165995.8, 165995.8, 298767.1
## $ geometry          <POINT [US_survey_foot]> POINT (1658494 182251.6), P…
## $ DEM               <dbl> 0.7807876, 1.1911737, 0.1216580

## mapping with leaflet
# https://rstudio.github.io/leaflet/
# https://cengel.github.io/R-spatial/mapping.html#web-mapping-with-leaflet


# endangered species occurrences in MD from GBIF
# shapefile at : https://github.com/collnell/GWU-visual/tree/master/MD_endangeredspecies
sps<-st_read('MD_endangeredspecies')

## Reading layer `MD_endangeredspecies' from data source `/Users/collnell/Dropbox/rstats/data viz/MD_endangeredspecies' using driver `ESRI Shapefile'
## Simple feature collection with 252 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 608948.3 ymin: 153422.9 xmax: 1867328 ymax: 754753.5
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=39.45 +lat_2=38.3 +lat_0=37.66666666666666 +lon_0=-77 +x_0=399999.9998983999 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs

glimpse(sps) #252 observations of various taxa

## Observations: 252
## Variables: 6
## $ species  <fct> Alasmidonta heterodon, Betula lenta, Betula lenta, Betu…
## $ kingdom  <fct> Animalia, Plantae, Plantae, Plantae, Plantae, Plantae, …
## $ phylum   <fct> Mollusca, Tracheophyta, Tracheophyta, Tracheophyta, Tra…
## $ order    <fct> Unionida, Fagales, Fagales, Fagales, Fagales, Fagales, …
## $ family   <fct> Unionidae, Betulaceae, Betulaceae, Betulaceae, Betulace…
## $ geometry <POINT [US_survey_foot]> POINT (1597053 438583.1), POINT (137…

st_crs(sps)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=lcc +lat_1=39.45 +lat_2=38.3 +lat_0=37.66666666666666 +lon_0=-77 +x_0=399999.9998983999 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

ggplot(sps)+
  geom_sf(aes(color=phylum), alpha=.5)

## see whether any endangered species have been observed near sites

# create buffer around sites
sites.buff<-st_buffer(sites.sf, dist=10000)
#buffering efffectively turns point datain to polygon

# plot buffered points
site.plot<-ggplot(dem.df)+
  geom_sf(data=co.bbox)+
  geom_sf(data=sites.buff, color='red')+
  theme_bw()

site.plot

# add points for species
sps.bbox<-sps[co.bbox,]
glimpse(sps.bbox)

## Observations: 5
## Variables: 6
## $ species  <fct> Pseudemys rubriventris, Pseudemys rubriventris, Pseudem…
## $ kingdom  <fct> Animalia, Animalia, Animalia, Animalia, Animalia
## $ phylum   <fct> Chordata, Chordata, Chordata, Chordata, Chordata
## $ order    <fct> Testudines, Testudines, Testudines, Testudines, Rodentia
## $ family   <fct> Emydidae, Emydidae, Emydidae, Emydidae, Sciuridae
## $ geometry <POINT [US_survey_foot]> POINT (1553412 301594.2), POINT (156…

site.plot+
  geom_sf(data=sps.bbox, color='blue')

# count points in polygons
st_within()
aggregate()

rnaturalearth

library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
glimpse(world)

## Observations: 241
## Variables: 64
## $ scalerank  <int> 3, 1, 1, 1, 1, 3, 3, 1, 1, 1, 3, 1, 5, 3, 1, 1, 1, 1,…
## $ featurecla <chr> "Admin-0 country", "Admin-0 country", "Admin-0 countr…
## $ labelrank  <dbl> 5, 3, 3, 6, 6, 6, 6, 4, 2, 6, 4, 4, 5, 6, 6, 2, 4, 5,…
## $ sovereignt <chr> "Netherlands", "Afghanistan", "Angola", "United Kingd…
## $ sov_a3     <chr> "NL1", "AFG", "AGO", "GB1", "ALB", "FI1", "AND", "ARE…
## $ adm0_dif   <dbl> 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0,…
## $ level      <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
## $ type       <chr> "Country", "Sovereign country", "Sovereign country", …
## $ admin      <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ adm0_a3    <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALD", "AND", "ARE…
## $ geou_dif   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ geounit    <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ gu_a3      <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALD", "AND", "ARE…
## $ su_dif     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ subunit    <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ su_a3      <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALD", "AND", "ARE…
## $ brk_diff   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ name       <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ name_long  <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ brk_a3     <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALD", "AND", "ARE…
## $ brk_name   <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ brk_group  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ abbrev     <chr> "Aruba", "Afg.", "Ang.", "Ang.", "Alb.", "Aland", "An…
## $ postal     <chr> "AW", "AF", "AO", "AI", "AL", "AI", "AND", "AE", "AR"…
## $ formal_en  <chr> "Aruba", "Islamic State of Afghanistan", "People's Re…
## $ formal_fr  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ note_adm0  <chr> "Neth.", NA, NA, "U.K.", NA, "Fin.", NA, NA, NA, NA, …
## $ note_brk   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Multiple…
## $ name_sort  <chr> "Aruba", "Afghanistan", "Angola", "Anguilla", "Albani…
## $ name_alt   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ mapcolor7  <dbl> 4, 5, 3, 6, 1, 4, 1, 2, 3, 3, 4, 4, 1, 7, 2, 1, 3, 1,…
## $ mapcolor8  <dbl> 2, 6, 2, 6, 4, 1, 4, 1, 1, 1, 5, 5, 2, 5, 2, 2, 1, 6,…
## $ mapcolor9  <dbl> 2, 8, 6, 6, 1, 4, 1, 3, 3, 2, 1, 1, 2, 9, 5, 2, 3, 5,…
## $ mapcolor13 <dbl> 9, 7, 1, 3, 6, 6, 8, 3, 13, 10, 1, NA, 7, 11, 5, 7, 4…
## $ pop_est    <dbl> 103065, 28400000, 12799293, 14436, 3639453, 27153, 83…
## $ gdp_md_est <dbl> 2258.0, 22270.0, 110300.0, 108.9, 21810.0, 1563.0, 36…
## $ pop_year   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ lastcensus <dbl> 2010, 1979, 1970, NA, 2001, NA, 1989, 2010, 2010, 200…
## $ gdp_year   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ economy    <chr> "6. Developing region", "7. Least developed region", …
## $ income_grp <chr> "2. High income: nonOECD", "5. Low income", "3. Upper…
## $ wikipedia  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ fips_10    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ iso_a2     <chr> "AW", "AF", "AO", "AI", "AL", "AX", "AD", "AE", "AR",…
## $ iso_a3     <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALA", "AND", "ARE…
## $ iso_n3     <chr> "533", "004", "024", "660", "008", "248", "020", "784…
## $ un_a3      <chr> "533", "004", "024", "660", "008", "248", "020", "784…
## $ wb_a2      <chr> "AW", "AF", "AO", NA, "AL", NA, "AD", "AE", "AR", "AM…
## $ wb_a3      <chr> "ABW", "AFG", "AGO", NA, "ALB", NA, "ADO", "ARE", "AR…
## $ woe_id     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ adm0_a3_is <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALA", "AND", "ARE…
## $ adm0_a3_us <chr> "ABW", "AFG", "AGO", "AIA", "ALB", "ALD", "AND", "ARE…
## $ adm0_a3_un <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ adm0_a3_wb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ continent  <chr> "North America", "Asia", "Africa", "North America", "…
## $ region_un  <chr> "Americas", "Asia", "Africa", "Americas", "Europe", "…
## $ subregion  <chr> "Caribbean", "Southern Asia", "Middle Africa", "Carib…
## $ region_wb  <chr> "Latin America & Caribbean", "South Asia", "Sub-Sahar…
## $ name_len   <dbl> 5, 11, 6, 8, 7, 5, 7, 20, 9, 7, 14, 10, 23, 22, 17, 9…
## $ long_len   <dbl> 5, 11, 6, 8, 7, 13, 7, 20, 9, 7, 14, 10, 27, 35, 19, …
## $ abbrev_len <dbl> 5, 4, 4, 4, 4, 5, 4, 6, 4, 4, 9, 4, 7, 10, 6, 4, 5, 4…
## $ tiny       <dbl> 4, NA, NA, NA, NA, 5, 5, NA, NA, NA, 3, NA, NA, 2, 4,…
## $ homepart   <dbl> NA, 1, 1, NA, 1, NA, 1, 1, 1, 1, NA, 1, NA, NA, 1, 1,…
## $ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((-69.89912 1..., MULTIPOL…

ggplot(data = world) +
  geom_sf()

ggplot(data = world) +
  geom_sf() +
  coord_sf(crs = "+init=epsg:3035")