3

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)?

Different roof types. My calculation so far only applies to pent- and saddle roofs.

  • 1
    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
2

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.

  • The project looks really promising. We have been working along the lines of something similar, but never got into the whole CityGML format since we started out with some LiDAR-redived DSM rasters. But CityGML might just be the right way to go, since the data availability is much better. – Achu Mani Mar 16 '17 at 10:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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