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

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)
# 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.
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)
}

``````

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 < O, 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)? • This is very interesting, but there are too many questions here -- there should, ideally, be just one. As far as Q2 is concerned, there are the ST_Slope and ST_Aspect functions in Postgis, but these would involve converting your GML to a raster, ie, it is not available on vectors, so there would be a significant preprocessing step there. Working on GML directly in databases is pretty painful in my experience -- Postgres, does support xpath, but I wouldn't want to write a processing pipeline using that. – John Powell Jul 12 '16 at 8:40
• The bonus question is very interesting, but you should probably post that as a separate question and tag it with geometry and the like. – John Powell Jul 12 '16 at 8:41
• @JohnBarça Do you know if it is possible to rasterize the CityGML in such a way that it outputs an elevation raster (ideally a DSM)? That would help me out a lot. For know, the only way I know how was to take the 3Dshape from the CityGML, convert it to TIN and the TIN then to raster (all in ArcGIS, which I would hate to rely on), with often faulty results. With regard to the Bonus question: I thought about it and now believe it should be solvable by calculating the normal vector of the polygon. Something along the line of this concept: stackoverflow.com/a/22838372/3189930 – Achu Mani Jul 12 '16 at 10:17
• Sorry, I haven't really done GML to raster conversion before . Going via TIN sounds like the best approach. Having to rely on ArcGIS is, indeed, unfortunate. – John Powell Jul 12 '16 at 13:02