1

I tried to do this over the weekend but had no luck - even after searching this forum. I'm trying to calculate the frontage (length of the property bordering the street) of a number of properties. I've got a shapefile of Hobart, Australia which I've slimmed down for this question as below.

Private Parcel = Houses

Casement = Streets

Example

#Download and unzip the file from https://drive.google.com/open?id=18yyrDSFz2jp0N3R0ehs80leVs9Qyv1wy
#libraries
library(tmap)
library(rgdal)
#Import
example <- readOGR(dsn="Example GIS",layer="Example")
#Visualise the data
tm_shape(example) + tm_polygons("CAD_TYPE1")
#I want to calculate the total length of each 'Private Parcel' polygon shared border with 'Casement'. I don't know how to do this.

Can anyone tell me how I might be able to do this?

After reading @Spacedman’s comment I noticed that some properties are aggregating.

#Previous Code
library(rgeos)
library(rgdal)
example <- readOGR(#File Location)
houses = example[example$CAD_TYPE1 == "Private Parcel",]
streets = example[example$CAD_TYPE1 == "Casement",]
frontages = gIntersection(houses, streets,byid=c(TRUE,FALSE))
length(frontages)
dim(houses)
plot(houses[12,])
plot(frontages[12],add=TRUE,lwd=4,col="red")
houses$frontage = gLength(frontages, byid=TRUE)

#My test
plot(houses)
plot(frontages,add=TRUE,lwd=4,col=houses@data$frontage)
#As you can see, the frontages are being aggregated between some properties.

enter image description here

#Here's my 2nd attempt
library(tidyverse)
library(sf)
example  = read_sf("Example GIS/Example.shp")
roads <- example[example$CAD_TYPE1=="Casement",]
house <- example[example$CAD_TYPE1=="Private Parcel",]
# intersection
int = st_intersection(house,roads)
# find out about the length of each line segment
int$len = st_length(int)
#Plot
ggplot(data=house) + geom_sf(fill=house$CID,alpha=0.5,col=NA) + geom_sf_label(aes(label=CID),size=2) + geom_sf(data=int,col=int$CID,size=2) 

enter image description here

Two properties are still aggregating - 1271016 and 1271014.

2
  • Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Please check out our short tour to learn about our focussed Q&A format. What happens when you run the code that you have presented?
    – PolyGeo
    Sep 9, 2019 at 2:36
  • 1
    Where did you look? All geometric operations on sp objects are handled by the rgeos package, and you can use gIntersection(house, roads) to get the intersection geometry - which for you should be the frontage. Loop over houses, and use gLength to get the line length. Some things can go wrong but try it and re-post when you get stuck.
    – Spacedman
    Sep 9, 2019 at 12:14

1 Answer 1

4

Split out the houses and the street features:

> houses = example[example$CAD_TYPE1 == "Private Parcel",]
> streets = example[example$CAD_TYPE1 == "Casement",]

Compute the intersection of houses and streets, using byid to compute the intersection for each house for all the street lines:

> frontages = gIntersection(houses, streets,byid=c(TRUE,FALSE))

That gives us the same number of frontage features as houses:

> length(frontages)
[1] 27
> dim(houses)
[1] 27 19

Lets see if this makes sense - pick a feature at random, number 12:

> plot(houses[12,])
> plot(frontages[12],add=TRUE,lwd=4,col="red")
> plot(streets, add=TRUE)

enter image description here

Looks good. If all you want is the length, extract this using gLength (again using byid to get the length of each of the 27 features) and add to the houses:

> houses$frontage = gLength(frontages, byid=TRUE)
> spplot(houses,"frontage")

enter image description here

Looks believable - the houses on the corners have more frontage.

Not sure if this code will have problems if, for example, the roads and houses don't exactly overlie, or if some house features have no frontage, or if houses cross street lines resulting in point intersections. But this works on your example and some adjustment might be needed once you scale up to more data...

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.