Calculating the distance between Ávila and Valencia
#Step1: Subsetting the citiesavila <-subset(esp1, NAME_2 =='Ávila')valencia <-subset(esp1, NAME_2 =='Valencia')#Step2: Union#This is necessary for later when we use st_combineavila <-st_union(avila)valencia <-st_union(valencia)#Step3: Calculating centroidsavila_ctr =st_centroid(avila)valencia_ctr =st_centroid(valencia)#Step4: Calculating distanced =st_distance(avila_ctr, valencia_ctr)d
st_buffer calculates buffers - polygons encompassing all points up to the specified distance from a given geometry.
We can for example calculate 250km buffers around Avila and Valencia.
To be precise about our distances, we will use EPSG codes
# Transform your data to a suitable projected CRS that uses meters as units# For example, you can use a UTM zone that covers your area of interest# Here, I'm assuming your data is around Avila, Spain, so I'm using UTM zone 30N (EPSG:2062)# See https://epsg.io/ and look for Spain# You should choose an appropriate UTM zone for your data. Note that otherwise, this would be in degrees#Step1: Using the relevant CRS to calculate distances in metersavila_ctr_sp <-st_transform(avila_ctr, crs =2062)valencia_ctr_sp <-st_transform(valencia_ctr, crs =2062)#Step2: Creating the buffer at 250km. Note that unit is in meters: st_crs(avila_ctr_sp) - LENGTHUNIT under Datumavila_250 =st_buffer(avila_ctr_sp, dist =250000)valencia_250 =st_buffer(valencia_ctr_sp, dist =250000)#Step3: Going back to the original crs in degreesavila_250 <-st_transform(avila_250, crs =st_crs(esp2))valencia_250 <-st_transform(valencia_250, crs =st_crs(esp2))
Buffers
And this is how we plot the shape file for the buffers that we have just created:
ggplot()+geom_sf(data=esp2)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=avila_250, color ="red", fill=NA)+geom_sf(data=valencia_250, color ="red", fill=NA)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Buffers
And this is what it looks like:
Geometry Difference from Pairs
There are four geometry-generating functions that operate on pairs of input geometries:
st_intersection
st_difference
st_sym_difference
st_union
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
union_area<-st_union(avila_250, valencia_250)ggplot()+geom_sf(data=esp2)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=union_area, color ="red", fill=NA)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
Geometry Difference from Pairs: union
Remember that st_union can also be applied to an individual layer,
This returns a dissolved union of all geometries:
Vector layer aggregation
Sometimes, we may not want to dissolve all features into a single geometry
Sometimes, we may want to dissolve by an attribute
Vector layer aggregation
Let us assume that we did not have the polygons for provinces in Spain.
We clearly do not need that many rows. We can aggregate all our railways
#Step2: Creating a unique variable to merge all the linesrails$all<-"all_rails"#Step3: Merging the linesrails2 <- rails %>%group_by(all) %>%summarise()
We now calculate the individual rail length per district
#Step5: Calculating individual rail length and turning it into kmintersection$indiv_rail<-st_length(intersection)intersection$indiv_rail<-set_units(intersection$indiv_rail, 'km')#Step6: Dividing rail length by polygon areaintersection$rail_density<-intersection$indiv_rail/intersection$area_sqkmintersection$rail_density<-as.numeric(intersection$rail_density)#Step7: Saving only relevant variablesintersection2<-subset(intersection, select=c("NAME_2", "rail_density"))#Step8: Dropping geometryintersection2<-st_drop_geometry(intersection2)#Step9: Performing a left mergeesp2c<-left_join(esp2b, intersection2, by =c("NAME_2"="NAME_2"))
Vector layer aggregation
We now calculate the individual rail length per district