CityGML: Calculate orientation and inclination angles from a polygon's LinearRing coordinates

Question

I have multiple CityGML files of different cities. Only one of them is in such a good condition that it includes the orientation and inclination angles of every roof polygon (as a gen:stringAttribute). Since this information is vital for my research, I want to extract the same information from the rest of the CityGML files.

With an XML parser for R, I should been able to extract every gml:LinearRing coordinate sequence that defines a roof polygon. I annotate the LinearRing coordinates inside a matrix and perform a calculation to determine the slope.

# A single bldg:RoofSurface's LinearRing (EPSG: 31467)
# (line breaks only added for readability in SE, revert 
# them to whitespaces when pasting this example into R)
gml.posList <- "4491058.587 5321312.714 524.130 
                4491066.325 5321319.284 524.130 
                4491064.062 5321322.034 521.847 
                4491056.270 5321315.460 521.827 
                4491058.587 5321312.714 524.130"

# We separate it at every whitespace and fill it into a three-columned matrix. 
mat.posList <- matrix(data = as.numeric(unlist(strsplit(gml.posList, split = " "))), 
ncol = 3, nrow = length(unlist(strsplit(gml.posList, split = " ")))/3,
byrow = T)

# The function starts a loop at the second row and calculates
# the distance (run) between each point and its neighboring 
# point in the row before. It then takes the elevation (rise)
# to calculate the slope. The slope-vector is then appended 
# to the input-matrix. 
addSlope <- function(mat.posList){
  slope <- NULL
  for(i in 2:nrow(mat.posList)){
    run <- sqrt(abs(mat.posList[i-1,1] - mat.posList[i,1])^2 + 
           (abs(mat.posList[i-1,2] - mat.posList[i,2])^2))
    rise <- abs(mat.posList[i-1,3] - mat.posList[i,3])
    slope[i] <- atan(rise/run) * 180 / pi
  }
  mat.posList <- cbind(mat.posList, slope)
  return(mat.posList)
}

mat.posList <- addSlope(mat.posList)

Now that I have determined the slope between the polygon vertices, I would love to to know the orientation / aspect angle at which the highest slope occurs, since this should (in most cases¹) be the orientation / aspect angle for the entire roof polygon. My first try, where I wanted top calculate the angle between an hypothetical north-vector and the vector between both points with the highest slope, works but is not precise.

# We determine the row at which the slope is the highest
j <- which(mat.posList[,4] == max(na.omit(mat.posList[,4])))

# we normalize the matrix which is necessary for the calculation in 
# the last line to work. For this, we remove the first six digits 
# of the coordinates, since they are the same for every vertice. 
# I know its messy, but I haven't come up with a more elegant solution.
normalized.mat <- matrix((as.numeric(substring(mat.posList[,1:2], 6)) - 
                  min(as.numeric(substring(mat.posList[,1:2], 6))))/ 
                  (max(as.numeric(substring(mat.posList[,1:2], 6))) - 
                  min(as.numeric(substring(mat.posList[,1:2], 6)))), ncol = 2)

# An O-rigin point, a point N-orth of it, and a point to which 
# we know the S-lope is the highest are created.
O <- unname(c(normalized.mat[j-1,1],  normalized.mat[j-1,2]))
N <- unname(c(normalized.mat[j-1,1], (normalized.mat[j-1,2]) + 
     (1 - normalized.mat[j-1,2])))
S <- unname(c(normalized.mat[j,1],  normalized.mat[j,2]))

# The angle between both vectors N-O and S-O is modeled after
# this response http://stackoverflow.com/a/1898026/3189930 . 
# The second line determines whether the angle falls inside 
# the 0 -180 or 180 to 360 range (whether S is east or west of O). 
orient <- acos(sum((N-O) * (S-O)) / ( sqrt(sum((N-O) * (N-O))) *
          sqrt(sum((S-O) * (S-O))))) * 180 / pi
orient <- ifelse(S[1] < O[1], 360 - orient, orient)

My questions:

1. How can I make the orientation angle more precise? The true value should be 319.7°. What I achieve is 39.45°, which, subtracted from 360° is close to the targeted orientation angle. But it is still of by almost an entire degree.
2. Can I do this directly in PostGIS/PostgreSQL? Since processing the entire CityGML Databases through R would take very long, and since I can't imagine that I am the first one who is interested in knowing the orientation and inclination angles of roof polygons, I would be very pleased to learn how to do this directly inside a CityGML database. I use 3DcityDB and its SPSHG-plugin to generate spreadsheets of the information stored inside the CityGML database, and would find it much easier to directly include the slope and orientation angles within those spreadsheets.

---

¹: When a roof is rectangular and the ridge and eave have the same length (a saddle roof). Which leads us to a bonus question: How should I perform such a calculation if a roof is hipped, or has an even more exotic shape (e.g., pyramid)?

Answer 1

This doesn't directly answer your questions, but some time ago I wrote a basic code to do what you're looking for. It extracts the orientation and inclination of RoofSurface polygons in CityGML, to estimate the solar irradiation of rooftops. The code is released on Github, so you might want to have a look. The part of the code that would probably interest you is here. The program creates a new enriched CityGML file with the information about each roof polygon. It's in Python and it's fairly simple, but I remember that it worked with large datasets.

Note that lots of polygons in publicly available CityGML data are not planar and contain other types of inconsistencies, hence the calculations may result in errors. My code has a built-in check for planarity to skip invalid polygons, but for more sophisticated validation you can use the free online validator val3dity.

By the way, with respect to your sentence:

Only one of them is in such a good condition that it includes the orientation and inclination angles of every roof polygon (as a gen:stringAttribute).

CityGML does not mandate attributes of the orientation of polygons and similar information, so the presence of such information doesn't say much about the condition of the data.