It is in the coordinate reference system of NAD83 and includes the full dataset for the US. You can map this in leaflet, but it is a big dataset and we want to focus on our spatial domains for projects, just the Inland Empire counties of Riverside and San Bernardino.
Thankfully, we learned how to apply a spatial filter using st_filter() last week as well.
IE_counties<-counties(state ='CA', cb =T, progress_bar =FALSE)|>#pull for California, lower detail level, year = 2022 ) filter(NAME%in%c('San Bernardino', 'Riverside'))
Simple feature collection with 2 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -117.8025 ymin: 33.42593 xmax: -114.1312 ymax: 35.80963
Geodetic CRS: NAD83
STATEFP COUNTYFP COUNTYNS AFFGEOID GEOID NAME
1 06 065 00277297 0500000US06065 06065 Riverside
2 06 071 00277300 0500000US06071 06071 San Bernardino
NAMELSAD STUSPS STATE_NAME LSAD ALAND AWATER
1 Riverside County CA California 06 18671846195 242807981
2 San Bernardino County CA California 06 51976311092 96418406
geometry
1 MULTIPOLYGON (((-117.6767 3...
2 MULTIPOLYGON (((-117.8025 3...
And now we have the two datasets we need to trim the rail data down to the local scale. Apply st_filter() to get an IE_rails dataset.
IE_rail<-st_filter(rails, IE_counties)head(IE_rail)# down to only 826 records from 120k
Simple feature collection with 6 features and 3 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: -117.6677 ymin: 33.73872 xmax: -116.2508 ymax: 34.00415
Geodetic CRS: NAD83
LINEARID FULLNAME MTFCC geometry
1 110458253217 Railroad R1011 LINESTRING (-116.2744 33.74...
2 110458253218 Railroad R1011 LINESTRING (-117.1677 34.00...
3 110458254923 At and Sf Rlwy R1011 LINESTRING (-117.5418 33.87...
4 110458254925 At and Sf Rlwy R1011 LINESTRING (-117.6677 33.87...
5 110458254934 At and Sf Rlwy R1011 LINESTRING (-117.2324 33.76...
6 110458254935 At and Sf Rlwy R1011 LINESTRING (-117.2334 33.76...
rm(ls =rails)# this removes the big rails dataset to free up memory, you can always reimport it.
Figure 23.1 shows the rail dataset for the Inland Empire. This is our basic map.
IE_rail|>st_transform(crs =4326)|># coordinate reference transformation to WGS84 for visualization leaflet()|>addTiles()|>setView(lat =34, lng =-117.3, zoom =10)|>addPolylines(color ='black', weight =4, dashArray ='2,5')|>addPolylines(data =IE_counties, color ='blue', weight =3)
Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'
Figure 23.1: Inland Empire rail links from tigris
Note
Adding fancy overlay options incoming!
Figure 23.2 shows the rail layer and a rail tile layer, and also displays the addLayersControl() function for allowing layers to be toggled on-off by the user.
IE_rail2<-IE_rail|>st_transform(crs =4326)leaflet()|>addTiles()|>#define group namesaddProviderTiles(providers$OpenRailwayMap, group ='RailTile')|>#this is a standard provider tileaddLayersControl(#This provides a toggle option for layers baseGroups =c('Basemap', 'RailTile'), #this is the names for tiles overlayGroups =c('RailTigris', 'County')#this is names for datalayers - polylines, polygons, markers, etc.)|>setView(lat =34, lng =-117.3, zoom =10)|>addPolylines(data =IE_rail2, color ='black', weight =4, dashArray ='2,5', group ='RailTigris')|>addPolylines(data =IE_counties, color ='blue', weight =3, group ='County')
Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
Need '+proj=longlat +datum=WGS84'
Figure 23.2: Inland Empire rail links with overlays!
23.3st_buffer()
The next function I want to demonstrate is the st_buffer() function. This allows you to create a user-defined distance from a geospatial object to help identify nearby features.
Let’s show an example using the IE_rail dataset.
IERailBuffer<-IE_rail|>st_buffer(dist =500, nQuadSegs =8)|># units are in the base units of the projection - could be meters or lat-lng degrees. NAD83 is a meter based projection, WGS84 is lat-lngst_simplify(dTolerance =150)|>#this smooths some wiggles and makes the data smaller bigger dTolerance would mean less wigglesst_transform(crs =4326)head(IERailBuffer)
Simple feature collection with 6 features and 3 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -117.6736 ymin: 33.73431 xmax: -116.2456 ymax: 34.0087
Geodetic CRS: WGS 84
LINEARID FULLNAME MTFCC geometry
1 110458253217 Railroad R1011 POLYGON ((-116.2746 33.7546...
2 110458253218 Railroad R1011 POLYGON ((-117.173 34.00299...
3 110458254923 At and Sf Rlwy R1011 POLYGON ((-117.5437 33.8720...
4 110458254925 At and Sf Rlwy R1011 POLYGON ((-117.6545 33.8744...
5 110458254934 At and Sf Rlwy R1011 POLYGON ((-117.2283 33.7684...
6 110458254935 At and Sf Rlwy R1011 POLYGON ((-117.2283 33.7684...
Note that IE_rail was a LINESTRING whereas IERailBuffer is now a POLYGON type.
Figure 23.3 shows a simple version of how it changed the data.