Archives de catégorie : 3D

[R] Buildings on Google Earth : from 2D to 3D



As I discover the marvels of (geo)R world, I found out this article (link)

It explains how to generate KML files from meuse, the source dataset delivered with sp package, using R software.

In the tutorial, KML files are created from scracth, writing each line of it.OGR driver only gives access to the options name and description (rdgal driver). So, « brute coding » may be necessary if you want to control other elements: color, heights for extruded objects…

Many people would like to visualize some buildings in 3D in google. The reference data that we have for buildings is often plain but there usually is a column specifying the height of each building. For example, for the french BDTOPO, the HAUTEUR column.

On the web, I found an arcScript for this but no opensource pre-built script or software.

Reading the article of the above mentioned tutorial called ploting 3D points gave me the idea to develop an R script for this purpose – which wasn’t so difficult, indeed.

Here is the complete script:

# import rgdal library which contains all necessary drivers
library(rgdal)
 
 
# shp import
 
buildings<-readOGR(dsn="C:/buildings.shp",layer="buildings")
 
 
#variable initialization
 
zColumn<-"HEIGHT"
 
 
filename<-"buildings.kml"
 
 
#header writing
 
write(paste("<?xml
 
version="1.0" encoding="ISO-8859-1" ?><kml
 
xmlns="http:#www.opengis.net/kml/2.2"><Document><Folder><name>buildingsge</name><Schema
 
name="buildingsge" id="buildingsge"></Schema>"
,
 
sep=""), filename)
 
 
#objects creation through a loop
 
for (x in 1:length(buildings@polygons)) {
 
 
write(paste("<Placemark><name>",
 
buildings@data$ID[x],
 
"</name><Style><LineStyle><color>7fbcbcbc</color></LineStyle><PolyStyle><color>7fbcbcbc</color></PolyStyle></Style><ExtendedData><SchemaData
 
schemaUrl="#buildingsge"><SimpleData name="nature">"
,
 
buildings@data$NATURE[x],"</SimpleData></SchemaData></ExtendedData><Polygon><extrude>1</extrude><altitudeMode>relativeToGround</altitudeMode><outerBoundaryIs><LinearRing><coordinates>",
 
sep=""), filename, append=TRUE)
 
 
# polygon coordinates injection
 
coordsPol<-paste(buildings@polygons[[x]]@Polygons[[1]]@coords[,1],
 
buildings@polygons[[x]]@Polygons[[1]]@coords[,2], buildings@data[[zColumn]][x],
 
sep=",")
 
write(coordsPol, filename, append=TRUE)
 
 
write("</coordinates></LinearRing></outerBoundaryIs></Polygon></Placemark>",filename,append=TRUE)}
 
 
# tail writing after exiting the loop
 
write("</Folder></Document></kml>",filename,append=TRUE)

Created by Pretty R at inside-R.org

Notes:
Here, I preferred to condense all the elements (<Placemark>, <coordinates>, etc…) in one line. Yet, it’s not optimized for reading. One may prefer individualizing each line by as many write() occurrences as elements like in the spatial analyst tuto.

writing

write("text",filename, append=TRUE)

is used to write the file from scratch

paste function
It’s used in order to concatenate several string elements 


Placemark name

buildings@data$ID[x]

each bulding is identified by its ID as extracted from the buildings shapefile. 

Coordinates

buildings@polygons[[x]]@Polygons[[1]]@coords[,1] 

is the X coordinate of the xth polygon

buildings@polygons[[x]]@Polygons[[1]]@coords[,2] 

is the Y

buildings@data[[zColumn]][x] 

accesses the zColumn value for the xth polygon/building.

Additional info

<ExtendedData<SchemaData
schemaUrl="#buildingsge">"nature">",
buildings@data$NATURE[x],"

will pop up the NATURE column value when clicking on an element (church, city hall, …)
If you’d like to understand the functions of some KML elementsn check the KML reference guide:
http://code.google.com/intl/fr/apis/kml/documentation/kmlreference.html

Importer le MNT de l’IGN avec QGIS & GRASS

Ici, nous allons voir comment importer le Modèle Numérique de Terrain de l’IGN dans QGIS/GRASS et comment visualiser le résultat dans le visionneur 3D NVIZ

Modification de l’en-tête du MNT d’origine
Le MNT de l’IGN se présente sous la forme d’un fichier à extension .asc comme ASCII

Si l’on ouvre le fichier dans un éditeur de texte, on trouve l’en-tête suivant qui donne la résolution, le nombre de colonnes/lignes composant, ainsi que les coordonnées X Y du coin bas gauche du raster:
ncols 3441
nrows 5401

xllcorner 1157987.50

yllcorner 6101987.50
cellsize 25.00

NODATA_value -9999

Pour être importé dans GRASS, le fichier .asc doit posséder la structure d’en-tête suivante:
cols: 3441
rows: 5401

west: 1157987.50
south: 6101987.50

north: 6237012.50

east: 1244012.50

Les éléments à calculer sont l’abcisse Xmax de l’extrêmité gauche (east) ainsi que l’ordonnée de l’extrêmité haute Ymax (north) du raster.

Le mode de calcul est très simple. Le voici:
east=ncols*cellsize+xllcorner
Ici:
east=(3441*25)+1157987.50
et pour north:
north=nrows*cellsize+yllcorner

En résumé, voici la correspondance entre l’ancien en-tête et le nouveau:
cols: ncols
rows: nrows

west: xllcorner

south: yllcorner

north: nrows*cellsize+yllcorner

east: ncols*cellsize+xllcorner

Une fois cela fait, enregistrez le MNT sous le nom
MNT_GRASS.asc par exemple

Import du fichier ASCII
Dans QGIS, ouvrez les outils GRASS et activez le module r.in.ascii.


Entrez les paramètres suivants:

Visualisation dans NVIZ
Une fois cela fait, vous pouvez visualiser votre MNT via le visionneur NVIZ mais avant le lancement de NVIZ, assurez-vous d’avoir bien défini votre région: l’étendue ainsi que la résolution de la zone d’intérêt à étudier..



Ici, nous avons utilisé la palette de couleur par défaut pour les rasters d’élévation. Il s’agit de la table de couleurs rainbow.

Drapage d’orthophotos, de scans
Je vous recommande fortement de lire le document suivant écrit par Martin Laloux qui traite des résolutions de raster dans NVIZ, dans les cas où l’on souhaite draper une orthophoto, un scan…
http://www.portailsig.org/content/grass-gis-region-resolution-resolution-3d-drapage-ou-l-apparente-quadrature-du-cercle

Ici, nous avons appliqué la méthode décrite ci-dessus en interpolant le MNT à la résolution 0.5 de l’orthophoto puis ayant appliqué la table de couleurs de l’orthophoto.



La valeur de fine a été fixée à 0.5, résolution de l’orthophoto.


Plus d’infos sur le traitement de rasters d’élévation, NVIZ:
r.in.ascii
r.resamp.interp
NVIZ