FrieseWoudloper/Geschiktheid panden provincie Groningen voor opwekken zonne-energie.RMD

## Geschiktheid panden provincie Groningen voor opwekken zonne-energie.RMD
---
title: "Welke provinciale panden zijn geschikt voor het opwekken van zonne-energie?"
author: "Willy Bakker"
date: "5 januari 2016"
output: html_document
---

## Vraagstelling en uitgangspunten
Een collega bij de Provincie Groningen vroeg zich af: Welke panden van de provincie Groningen zijn geschikt voor het opwekken van zonne-energie? Een interessante vraag! Bovendien een uitgelezen kans om me verder te verdiepen in geografische analyses met R. In de kerstvakantie ben ik er mee aan de slag gegaan. Dit document bevat de code, documentatie en resultaten van mijn analyse.

Mijn uitgangspunt was dat de analyse van begin tot eind reproduceerbaar zou moeten zijn. Daaruit voortvloeiend heb ik de volgende regels gehanteerd:

* Maak gebruik van open data.
* Maak gebruik van open source software (R).
* Genereer zowel de code als de documentatie vanuit ��n Markdown-document.
* Deel de alles via het internet (GitHub).

Voor mijn analyse heb ik de volgende databronnen gebruikt:

* [Gebouwen](http://www.provinciegroningen.nl/objecten/gebouwen/) van de provincie Groningen
* [Zonatlas](http://www.zonatlas.nl/home/)
* [Basisregistratie Adressen en Gebouwen](https://www.pdok.nl/nl/service/wfs-bag) (BAG)
* [Bestuurlijke grenzen](https://www.pdok.nl/nl/service/wfs-bestuurlijke-grenzen)

## Analyse

### Voorbereiding: Importeren packages
Om te beginnen importeer ik de packages die nodig zijn voor de analyse. Dat is een flinke lijst ;-)
```{r, messages = FALSE, warning = FALSE, results = 'hide'}
library(httr)
library(xlsx)
library(sp)
library(rgdal)
library(rgeos)
library(maptools)
library(png)
library(raster)
library(plyr)
```

### Stap 1: Inlezen BAG-identificatienummers van panden en verblijfsobjecten
Vervolgens download ik het Excel-bestand Gebouwen uit het objectenportaal van de provincie Groningen. Voor de analyse ben ik alleen ge�nteresseerd in de kolommen Complex en BAGID.

```{r, messages = FALSE, warning = FALSE, results = 'hide'}
url <- "http://geoservices.provinciegroningen.nl/bestanden/ObjectenBlauw/gebouwen/DATA/gebouwen.xls"
tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".xls")
download.file(url, destfile = tmp, mode = "wb")
rijen <- read.xlsx2(tmp, sheetIndex = 1, colIndex = c(3, 6), colClasses = "character", stringsAsFactors = FALSE)
identificatienrs <- unique(rijen[rijen$BAGID != "", c("BAGID")])
unlink(tmp)
```

In het Excel-bestand staan panden, verblijfsobjecten en overige objecten door elkaar. Helaas kan ik de laatste categorie niet meenemen in de analyse, omdat de BAG voor deze objecten geen geometrie bevat. Er is voor deze gegevens ook geen alternatieve bron beschikbaar.

Sommige panden of verblijfsobjecten komen meerdere malen voor in het Excel-bestand. Deze dubbele waarden filter ik uit de dataset.

Uiteindelijk neem ik `r length(identificatienrs)` van de in totaal `r nrow(rijen)` rijen in het Excel-bestand mee in de analyse.

### Stap 2: Omzetten verblijfsobjecten naar panden
In de lijst met BAG-identificatienummers staan panden en verblijfsobjecten door elkaar heen. Dat is lastig, omdat verblijfsobjecten in de BAG meestal een puntgeometrie hebben. Panden zijn vlakken. Ik heb vlakken nodig om te bepalen of het dak van een gebouw geschikt is voor zonne-energie. Ik moet dus eerst de verblijfsobjecten in de lijst omzetten naar panden.
Met behulp van de BAG kan ik de identificatienummers opvragen van panden die horen bij een specifiek verblijfsobject. Hierbij gelden de volgende aandachtspunten:

* E�n verblijfsobject kan gekoppeld zijn aan meerdere panden.
* Als een pand meerdere gebruiksdoelen heeft, retourneert de BAG-service hetzelfde pand meerdere malen: ��n keer voor ieder gebruiksdoel.

```{r}
base.url <- "http://geodata.nationaalgeoregister.nl/bag/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bag:verblijfsobject&outputFormat=json&srsName=epsg:4326"

identificatienrs.panden <- sapply(identificatienrs, function(x){
    url <- paste0(base.url, "&cql_filter=(bag:identificatie='", x, "')")
    tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
    id <- character()
    r <- content(GET(url, write_disk(tmp, overwrite = TRUE)), 'parse')
    if (r$totalFeatures > 0) {
        # Het object is een verblijfsobject. Lees de bijbehorende pand-ID's in uit de response van de BAG-service.
        for (i in 1:r$totalFeatures) id <- c(id, r$features[[i]]$properties$pandidentificatie)
    }
    else {
        # Het object is een pand. Geen actie nodig. Bewaar het pand-ID uit de BAG.
        id <- x
    }
    unlink(tmp)
    return (id)
})

identificatienrs.panden <- unique(unlist(identificatienrs.panden, use.names = FALSE))
```

### Stap 3: Inlezen gemeentegrens Groningen
Er zijn twee versies van de Zonatlas: ��n voor de gemeente Groningen en ��n voor de rest van de provincie. Om voor een pand de juiste gegevens uit de Zonatlas te selecteren, moet ik eerst weten of het pand binnen of buiten de gemeente Groningen ligt. Daarvoor heb ik de gemeentegrens nodig.

In onderstaande code vraag ik expliciet aan de BAG-service om een GeoJSON-bestand retour te sturen. Het GeoJSON-bestand is een workaround. Het lukt me namelijk niet om de reguliere GML-respons van de WFS-service te verwerken in R. Het is een probleem waar ook [andere Windows-gebruikers](http://stackoverflow.com/questions/32034495/list-aviable-wfs-layers-and-read-into-data-frame-with-rgdal) tegenaan lopen.

Om fouten verderop in de analyse te voorkomen, zorg ik ervoor dat alle geografische databronnen worden getransformeerd naar [EPSG:4326](http://spatialreference.org/ref/epsg/wgs-84/). R geeft een waarschuwing wanneer ik de gemeentegrenzen omzet van WGS84 naar EPSG:4326. Deze melding negeer ik. In theorie zijn WGS84 en EPSG:4326 namelijk gelijk aan elkaar, maar voor de zekerheid voer ik de transformatie toch maar uit.

```{r, message = FALSE, results = 'hide'}
url <- "http://geodata.nationaalgeoregister.nl/bestuurlijkegrenzen/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bestuurlijkegrenzen:gemeenten&outputFormat=json&srsName=epsg:4326&cql_filter=(gemeentenaam=%27Groningen%27)"
tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
r <- GET(url, write_disk(tmp, overwrite = TRUE))
gemeentegrens.Groningen <- readOGR(tmp, "OGRGeoJSON")
proj4string(gemeentegrens.Groningen) <- CRS("+init=epsg:4326")
unlink(tmp)
```

### Stap 4: Inlezen geometrie van panden
De volgende stap is om aan de lijst met pandidentificatienummers een kolom met de bijbehorende geometrie toe te voegen. Ik vraag per pand de geometrie op bij de BAG WFS-service. Hierbij moet ik opnieuw met een aantal zaken rekening houden:

* Tranformeer de geometrie naar EPSG:4326.
# Negeer de waarschuwingen "Z-dimension discarded" en "A new CRS was assigned to an object with an existing CRS".
* Filter dubbele rijen (die ontstaan zijn door meerdere gebruikersdoelen voor ��n pand) uit het resultaat.
* Zorg er voor dat iedere rij ('feature') in het resultaat een uniek ID heeft. Dit ID is gelijk aan het pandidentificatienummer.
* Let op dat het pandidentificatienummer correct is! De service retourneert het nummer niet als karakterstring, maar als numerieke waarde. Hierdoor gaan voorloopnullen verloren.
* Voeg een boolean attribuut toe waarvan de waarde aangeeft of de betreffende feature binnen of buiten de gemeente Groningen ligt.


Eerst maak ik voor ieder pand een SpatialPolygonsDataFrame (SPDF) aan, op basis van de gegevens die de BAG-service voor het pand teruggeeft.
```{r, message = FALSE, warning = FALSE, results = 'hide'}
base.url <- "http://geodata.nationaalgeoregister.nl/bag/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bag:pand&outputFormat=json&srsName=epsg:4326"
panden <- lapply(identificatienrs.panden, function(x){
    url <- paste0(base.url, "&cql_filter=(bag:identificatie=", "'", x, "')")
    tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
    r <- content(GET(url, write_disk(tmp, overwrite = TRUE)), 'parse')
    spdf <- readOGR(tmp, "OGRGeoJSON", stringsAsFactors = FALSE)
    proj4string(spdf) <- CRS("+init=epsg:4326")
    spdf <- spdf[!duplicated(spdf$identificatie),]
    id <- paste0(c(rep("0", 16 - nchar(spdf$identificatie)), spdf$identificatie), collapse = "")
    spdf <- spChFIDs(spdf, id)
    spdf$pand.ID <- id
    spdf$in.gemeente.Groningen <- !all(is.na(over(gemeentegrens.Groningen, spdf)))
    unlink(tmp)
    return (spdf[, c("pand.ID", "in.gemeente.Groningen")])
})
```
Vervolgens voeg ik de lijst met SPDF's samen tot ��n groot SPDF met meerdere panden.
```{r}
panden <- Reduce(spRbind, panden)
```

De resulterende SpatialPolygonsDataFrame bevat voor ieder pand nu naast de geometrie ook twee attributen: het pandidentificatienummer uit de BAG en een boolean die aangeeft of het pand wel of niet in de gemeente Groningen ligt.
De rijnamen zijn gelijk aan de pandidentificatienummers uit de BAG.
```{r, echo = FALSE}
head(panden)
```

### Stap 5: Bepaal welke tegels uit de Zonatlas nodig zijn voor de analyse
De gegevens in de Zonatlas komen van een Tile Map Service (TMS). Het tiling schema van de Zonatlas TMS is gelijk aan dat van OpenStreetMap, alleen zijn de y-waarden van de tegels ('tiles') anders genummerd. Bovendien zijn ze af- in plaats van oplopend.

In de [OpenStreetMap Wiki](http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames) is beschreven hoe je op basis van een co�rdinaat en het gewenste zoomniveau de URL van de bijbehorende OpenStreetMap tegel kunt bepalen. De code uit de OpenStreetMap wiki heb ik aangepast, zodat onderstaande functie niet de x- en y- waarden van de OpenStreetMap tegel, maar die van de Zonatlas tegel retourneert.

```{r}
deg2num <- function(lat, lon, zoom) {
    x <- floor((lon + 180.0) / 360.0 * 2.0 ^ zoom)
    y = 262143 - floor((1.0 - log(tan(lat * pi / 180) + (1 / cos(lat * pi / 180))) / pi) / 2.0 * 2.0 ^ zoom) # LET OP: correctie t.o.v. OSM tile schema!
    return (c(x, y))
}
```

In deze stap maak ik twee lijsten:

* Een ontdubbelde lijst met de URL's van alle tegels die gedownload moeten worden voor de analyse. (Strikt genomen is dit een data frame, en geen lijst.)
* Een lijst met alle tegels per pand. In deze lijst zitten wel dubbelingen.

Let op: het eerste deel van de URL is voor panden binnen en buiten de gemeente Groningen verschillend!

```{r, message = FALSE, warning = FALSE, results = 'hide'}
zoom <- 18
tiles <- data.frame()
tiles.per.pand <- list()

for (k in 1:length(panden)){
    bb <- bbox(panden@polygons[[k]])
    lb <- deg2num(bb["y","max"], bb["x","min"], zoom) # linksboven
    ro <- deg2num(bb["y","min"], bb["x","max"], zoom) # rechtsonder
    l <- list()
    for (y in lb[2]:ro[2]){
        for (x in ro[1]:lb[1]){
            if (panden$in.gemeente.Groningen[[k]]) base.url <- "http://viewermaia.zonatlas.nl/groningen/layers/pv/"
            else base.url <- "http://viewermaia.zonatlas.nl/prov_groningen/layers/pv/"
            l[[length(l)+ 1]] <- paste0(x, "-", y, ".grd")
            tiles <- rbind(tiles, data.frame(x, y, paste0(base.url, zoom, "/", x,"/", y, ".png"), stringsAsFactors = FALSE))
        }
    }
    tiles.per.pand[[length(tiles.per.pand) + 1]] <- c(panden$pand.ID[[k]], list(l))
}

tiles <- unique(tiles)
colnames(tiles) = c("x", "y", "url")
```

### Stap 5: Download de Zonatlas tegels en zet kleuren om in waarden
Voordat je de Zonatlas tegels kunt downloaden en opslaan voor verdere analyse, moet je ze eerst georefereren. Hiervoor gebruik ik een aantal functies afkomstig uit de OpenStreetMap Wiki. Waar nodig heb ik de code aangepast op het tiling schema van de Zonatlas.

```{r}
tile2lon <- function(x, zoom) {
    return (x / 2.0 ^ zoom * 360.0 - 180.0)
}

tile2lat <- function(y, zoom) {
    y.OSM <- 262143 - y # LET OP: correctie t.o.v. OSM tile schema!
    n <- 2.0 ^ zoom
    lat_rad <- atan(sinh(pi * (1 - 2 * y.OSM / n)))
    lat_deg <- 180.0 * (lat_rad / pi)
    return (lat_deg)
}

get.bounding.box <- function(x, y, zoom){
    return (c(ymin = tile2lat(y - 1, zoom), ymax = tile2lat(y, zoom), xmin = tile2lon(x, zoom), xmax = tile2lon(x + 1, zoom)))
}
```

Nu kunnen de Zonatlas tegels gedownload worden. Het zijn PNG-bestanden, die worden ingelezen en opgeslagen in het native raster formaat. Dit geeft me de mogelijkheid om kleuren in het rasterbestand om te zetten naar waarden (rood = minder geschikt, geel = geschikt, groen = zeer geschikt).

```{r, message = FALSE, warning = FALSE, results = 'hide'}
for (i in 1:nrow(tiles)) {
    url <- tiles$url[i]
    tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".png")
    download.file(url, destfile = tmp, mode = "wb")

    if (file.size(tmp) > 222) {
        img <- as.numeric(readPNG(tmp, native = TRUE))
        img[img %in% c(-14737455, -14737456)] <- 1 # rood  = minder geschikt
        img[img == -14423829] <- 2                 # geel  = geschikt
        img[img == -14114800] <- 3                 # groen = zeer geschikt
        img[!img %in% c(1, 2, 3)] <- NA
        bbox <- get.bounding.box(tiles$x[i], tiles$y[i], zoom)
        r <- raster(matrix(img, nrow = 256, ncol = 256, byrow = TRUE),
                    xmn = bbox["xmin"], xmx = bbox["xmax"], ymn = bbox["ymin"], ymx = bbox["ymax"],
                    crs = "+init=epsg:4326")
        writeRaster(r, paste0(tiles$x[i], "-", tiles$y[i], ".grd"), format = "raster", overwrite = TRUE)
    }
    unlink(tmp)
}
```

### Stap 6: Bepaal voor ieder pand de waarde in de Zonatlas
In de laatste stap voer ik per pand de volgende deelstappen uit:
1. Als er meerdere Zonatlas tegels bij het pand horen, voeg de native rasterbestanden dan samen tot ��n.
2. Selecteer alleen de rasterwaarden die binnen de begrenzing van het pand vallen.
3. Bereken de mediaan van de rasterwaarden binnen de begrenzing van het pand.
4. Decodeer de mediaan naar een waarde voor de geschiktheid van het pand voor het opwekken van zonne-energie

Let op: de waarde NA geeft aan dat de geschiktheid van het pand voor het opwekken van zonne-energie niet bepaald kon worden.

```{r, message = FALSE, warning = FALSE, results = 'hide'}
decode.color <- function (x){
    if (is.na(x)) return (NA)
    else if (x == 1) return ("minder geschikt")
    else if (x == 2) return ("geschikt")
    else if (x == 3) return ("zeer geschikt")
}

for (i in 1:length(tiles.per.pand)) {
    tiles <- tiles.per.pand[[i]][[2]]
    tiles <- tiles[sapply(tiles, file.exists)]
    if (length(tiles) == 0) {
        # Voor deze panden bevat de Zonatlas geen informatie.
        panden[tiles.per.pand[[i]][[1]], c("zonenerg")] <- NA
    } else {
        if (length(tiles) == 1) {
            r <- raster(tiles[[1]])
        } else if (length(tiles) > 1) {
            r <- do.call(merge, lapply(tiles, function(x){t <- raster(x); origin(t) <- c(0,0); return (t) }))
        }
        sp <- panden[tiles.per.pand[[i]][[1]],]
        panden[tiles.per.pand[[i]][[1]], c("zonenerg")] <- decode.color(summary(mask(r, sp))[3])
    }
}
```

Het eindresultaat schrijf ik tenslotte weg naar een shape-bestand met de naam panden-provincie.
In een shape-bestand mogen attribuutnamen niet langer zijn dan 8 karakters, vandaar de lelijke namen.
```{r, message = FALSE, warning = FALSE, results = 'hide'}
panden <- spTransform(panden[, c("pand.ID",  "zonenerg")], CRS("+init=epsg:28992"))
sapply(c("panden-provincie.dbf", "panden-provincie.prj", "panden-provincie.shp", "panden-provincie.shx"), unlink)
writeOGR(panden, dsn = getwd(), layer ='panden-provincie', driver = 'ESRI Shapefile')
```

## Resultaten
Van de `r nrow(rijen)` objecten in de spreadsheet Gebouwen.xls uit het projecten portaal, zijn `r length(identificatienrs)` meegenomen in de analyse.
Voor `r (nrow(rijen) - length(identificatienrs))` kon geen vlakgeometrie bepaald worden, vandaar dat ze niet meedoen.

kable(head(panden), caption = "De eerste 6 rijen", row.names = FALSE, col.names = c("Pand ID", "Geschiktheid voor opwekken zonne-energie"), padding = 2)
table(panden$zonenerg, useNA = "ifany")

panden <- readOGR(dsn = getwd(), layer ='panden-provincie')
for (i in 1:nrow(panden)) {
        panden$oppervlakte.m2[i] <- panden@polygons[[i]]@area
}

x <- arrange(panden@data[panden$zonenerg == 'zeer geschikt',], desc(oppervlakte.m2))
x <- rename(x, c("pand_ID" = "BAGID"))
data <- join(x, rijen, by = "BAGID")
head(data[, c("BAGID", "Complex", "zonenerg", "oppervlakte.m2")], n = 5)
	---
	title: "Welke provinciale panden zijn geschikt voor het opwekken van zonne-energie?"
	author: "Willy Bakker"
	date: "5 januari 2016"
	output: html_document
	---

	## Vraagstelling en uitgangspunten
	Een collega bij de Provincie Groningen vroeg zich af: Welke panden van de provincie Groningen zijn geschikt voor het opwekken van zonne-energie? Een interessante vraag! Bovendien een uitgelezen kans om me verder te verdiepen in geografische analyses met R. In de kerstvakantie ben ik er mee aan de slag gegaan. Dit document bevat de code, documentatie en resultaten van mijn analyse.

	Mijn uitgangspunt was dat de analyse van begin tot eind reproduceerbaar zou moeten zijn. Daaruit voortvloeiend heb ik de volgende regels gehanteerd:

	* Maak gebruik van open data.
	* Maak gebruik van open source software (R).
	* Genereer zowel de code als de documentatie vanuit ��n Markdown-document.
	* Deel de alles via het internet (GitHub).

	Voor mijn analyse heb ik de volgende databronnen gebruikt:

	* [Gebouwen](http://www.provinciegroningen.nl/objecten/gebouwen/) van de provincie Groningen
	* [Zonatlas](http://www.zonatlas.nl/home/)
	* [Basisregistratie Adressen en Gebouwen](https://www.pdok.nl/nl/service/wfs-bag) (BAG)
	* [Bestuurlijke grenzen](https://www.pdok.nl/nl/service/wfs-bestuurlijke-grenzen)

	## Analyse

	### Voorbereiding: Importeren packages
	Om te beginnen importeer ik de packages die nodig zijn voor de analyse. Dat is een flinke lijst ;-)
	```{r, messages = FALSE, warning = FALSE, results = 'hide'}
	library(httr)
	library(xlsx)
	library(sp)
	library(rgdal)
	library(rgeos)
	library(maptools)
	library(png)
	library(raster)
	library(plyr)
	```

	### Stap 1: Inlezen BAG-identificatienummers van panden en verblijfsobjecten
	Vervolgens download ik het Excel-bestand Gebouwen uit het objectenportaal van de provincie Groningen. Voor de analyse ben ik alleen ge�nteresseerd in de kolommen Complex en BAGID.

	```{r, messages = FALSE, warning = FALSE, results = 'hide'}
	url <- "http://geoservices.provinciegroningen.nl/bestanden/ObjectenBlauw/gebouwen/DATA/gebouwen.xls"
	tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".xls")
	download.file(url, destfile = tmp, mode = "wb")
	rijen <- read.xlsx2(tmp, sheetIndex = 1, colIndex = c(3, 6), colClasses = "character", stringsAsFactors = FALSE)
	identificatienrs <- unique(rijen[rijen$BAGID != "", c("BAGID")])
	unlink(tmp)
	```

	In het Excel-bestand staan panden, verblijfsobjecten en overige objecten door elkaar. Helaas kan ik de laatste categorie niet meenemen in de analyse, omdat de BAG voor deze objecten geen geometrie bevat. Er is voor deze gegevens ook geen alternatieve bron beschikbaar.

	Sommige panden of verblijfsobjecten komen meerdere malen voor in het Excel-bestand. Deze dubbele waarden filter ik uit de dataset.

	Uiteindelijk neem ik `r length(identificatienrs)` van de in totaal `r nrow(rijen)` rijen in het Excel-bestand mee in de analyse.

	### Stap 2: Omzetten verblijfsobjecten naar panden
	In de lijst met BAG-identificatienummers staan panden en verblijfsobjecten door elkaar heen. Dat is lastig, omdat verblijfsobjecten in de BAG meestal een puntgeometrie hebben. Panden zijn vlakken. Ik heb vlakken nodig om te bepalen of het dak van een gebouw geschikt is voor zonne-energie. Ik moet dus eerst de verblijfsobjecten in de lijst omzetten naar panden.
	Met behulp van de BAG kan ik de identificatienummers opvragen van panden die horen bij een specifiek verblijfsobject. Hierbij gelden de volgende aandachtspunten:

	* E�n verblijfsobject kan gekoppeld zijn aan meerdere panden.
	* Als een pand meerdere gebruiksdoelen heeft, retourneert de BAG-service hetzelfde pand meerdere malen: ��n keer voor ieder gebruiksdoel.

	```{r}
	base.url <- "http://geodata.nationaalgeoregister.nl/bag/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bag:verblijfsobject&outputFormat=json&srsName=epsg:4326"

	identificatienrs.panden <- sapply(identificatienrs, function(x){
	url <- paste0(base.url, "&cql_filter=(bag:identificatie='", x, "')")
	tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
	id <- character()
	r <- content(GET(url, write_disk(tmp, overwrite = TRUE)), 'parse')
	if (r$totalFeatures > 0) {
	# Het object is een verblijfsobject. Lees de bijbehorende pand-ID's in uit de response van de BAG-service.
	for (i in 1:r$totalFeatures) id <- c(id, r$features[[i]]$properties$pandidentificatie)
	}
	else {
	# Het object is een pand. Geen actie nodig. Bewaar het pand-ID uit de BAG.
	id <- x
	}
	unlink(tmp)
	return (id)
	})

	identificatienrs.panden <- unique(unlist(identificatienrs.panden, use.names = FALSE))
	```

	### Stap 3: Inlezen gemeentegrens Groningen
	Er zijn twee versies van de Zonatlas: ��n voor de gemeente Groningen en ��n voor de rest van de provincie. Om voor een pand de juiste gegevens uit de Zonatlas te selecteren, moet ik eerst weten of het pand binnen of buiten de gemeente Groningen ligt. Daarvoor heb ik de gemeentegrens nodig.

	In onderstaande code vraag ik expliciet aan de BAG-service om een GeoJSON-bestand retour te sturen. Het GeoJSON-bestand is een workaround. Het lukt me namelijk niet om de reguliere GML-respons van de WFS-service te verwerken in R. Het is een probleem waar ook [andere Windows-gebruikers](http://stackoverflow.com/questions/32034495/list-aviable-wfs-layers-and-read-into-data-frame-with-rgdal) tegenaan lopen.

	Om fouten verderop in de analyse te voorkomen, zorg ik ervoor dat alle geografische databronnen worden getransformeerd naar [EPSG:4326](http://spatialreference.org/ref/epsg/wgs-84/). R geeft een waarschuwing wanneer ik de gemeentegrenzen omzet van WGS84 naar EPSG:4326. Deze melding negeer ik. In theorie zijn WGS84 en EPSG:4326 namelijk gelijk aan elkaar, maar voor de zekerheid voer ik de transformatie toch maar uit.

	```{r, message = FALSE, results = 'hide'}
	url <- "http://geodata.nationaalgeoregister.nl/bestuurlijkegrenzen/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bestuurlijkegrenzen:gemeenten&outputFormat=json&srsName=epsg:4326&cql_filter=(gemeentenaam=%27Groningen%27)"
	tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
	r <- GET(url, write_disk(tmp, overwrite = TRUE))
	gemeentegrens.Groningen <- readOGR(tmp, "OGRGeoJSON")
	proj4string(gemeentegrens.Groningen) <- CRS("+init=epsg:4326")
	unlink(tmp)
	```

	### Stap 4: Inlezen geometrie van panden
	De volgende stap is om aan de lijst met pandidentificatienummers een kolom met de bijbehorende geometrie toe te voegen. Ik vraag per pand de geometrie op bij de BAG WFS-service. Hierbij moet ik opnieuw met een aantal zaken rekening houden:

	* Tranformeer de geometrie naar EPSG:4326.
	# Negeer de waarschuwingen "Z-dimension discarded" en "A new CRS was assigned to an object with an existing CRS".
	* Filter dubbele rijen (die ontstaan zijn door meerdere gebruikersdoelen voor ��n pand) uit het resultaat.
	* Zorg er voor dat iedere rij ('feature') in het resultaat een uniek ID heeft. Dit ID is gelijk aan het pandidentificatienummer.
	* Let op dat het pandidentificatienummer correct is! De service retourneert het nummer niet als karakterstring, maar als numerieke waarde. Hierdoor gaan voorloopnullen verloren.
	* Voeg een boolean attribuut toe waarvan de waarde aangeeft of de betreffende feature binnen of buiten de gemeente Groningen ligt.


	Eerst maak ik voor ieder pand een SpatialPolygonsDataFrame (SPDF) aan, op basis van de gegevens die de BAG-service voor het pand teruggeeft.
	```{r, message = FALSE, warning = FALSE, results = 'hide'}
	base.url <- "http://geodata.nationaalgeoregister.nl/bag/wfs?service=wfs&version=2.0.0&request=GetFeature&typeName=bag:pand&outputFormat=json&srsName=epsg:4326"
	panden <- lapply(identificatienrs.panden, function(x){
	url <- paste0(base.url, "&cql_filter=(bag:identificatie=", "'", x, "')")
	tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".json")
	r <- content(GET(url, write_disk(tmp, overwrite = TRUE)), 'parse')
	spdf <- readOGR(tmp, "OGRGeoJSON", stringsAsFactors = FALSE)
	proj4string(spdf) <- CRS("+init=epsg:4326")
	spdf <- spdf[!duplicated(spdf$identificatie),]
	id <- paste0(c(rep("0", 16 - nchar(spdf$identificatie)), spdf$identificatie), collapse = "")
	spdf <- spChFIDs(spdf, id)
	spdf$pand.ID <- id
	spdf$in.gemeente.Groningen <- !all(is.na(over(gemeentegrens.Groningen, spdf)))
	unlink(tmp)
	return (spdf[, c("pand.ID", "in.gemeente.Groningen")])
	})
	```
	Vervolgens voeg ik de lijst met SPDF's samen tot ��n groot SPDF met meerdere panden.
	```{r}
	panden <- Reduce(spRbind, panden)
	```

	De resulterende SpatialPolygonsDataFrame bevat voor ieder pand nu naast de geometrie ook twee attributen: het pandidentificatienummer uit de BAG en een boolean die aangeeft of het pand wel of niet in de gemeente Groningen ligt.
	De rijnamen zijn gelijk aan de pandidentificatienummers uit de BAG.
	```{r, echo = FALSE}
	head(panden)
	```

	### Stap 5: Bepaal welke tegels uit de Zonatlas nodig zijn voor de analyse
	De gegevens in de Zonatlas komen van een Tile Map Service (TMS). Het tiling schema van de Zonatlas TMS is gelijk aan dat van OpenStreetMap, alleen zijn de y-waarden van de tegels ('tiles') anders genummerd. Bovendien zijn ze af- in plaats van oplopend.

	In de [OpenStreetMap Wiki](http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames) is beschreven hoe je op basis van een co�rdinaat en het gewenste zoomniveau de URL van de bijbehorende OpenStreetMap tegel kunt bepalen. De code uit de OpenStreetMap wiki heb ik aangepast, zodat onderstaande functie niet de x- en y- waarden van de OpenStreetMap tegel, maar die van de Zonatlas tegel retourneert.

	```{r}
	deg2num <- function(lat, lon, zoom) {
	x <- floor((lon + 180.0) / 360.0 * 2.0 ^ zoom)
	y = 262143 - floor((1.0 - log(tan(lat * pi / 180) + (1 / cos(lat * pi / 180))) / pi) / 2.0 * 2.0 ^ zoom) # LET OP: correctie t.o.v. OSM tile schema!
	return (c(x, y))
	}
	```

	In deze stap maak ik twee lijsten:

	* Een ontdubbelde lijst met de URL's van alle tegels die gedownload moeten worden voor de analyse. (Strikt genomen is dit een data frame, en geen lijst.)
	* Een lijst met alle tegels per pand. In deze lijst zitten wel dubbelingen.

	Let op: het eerste deel van de URL is voor panden binnen en buiten de gemeente Groningen verschillend!

	```{r, message = FALSE, warning = FALSE, results = 'hide'}
	zoom <- 18
	tiles <- data.frame()
	tiles.per.pand <- list()

	for (k in 1:length(panden)){
	bb <- bbox(panden@polygons[[k]])
	lb <- deg2num(bb["y","max"], bb["x","min"], zoom) # linksboven
	ro <- deg2num(bb["y","min"], bb["x","max"], zoom) # rechtsonder
	l <- list()
	for (y in lb[2]:ro[2]){
	for (x in ro[1]:lb[1]){
	if (panden$in.gemeente.Groningen[[k]]) base.url <- "http://viewermaia.zonatlas.nl/groningen/layers/pv/"
	else base.url <- "http://viewermaia.zonatlas.nl/prov_groningen/layers/pv/"
	l[[length(l)+ 1]] <- paste0(x, "-", y, ".grd")
	tiles <- rbind(tiles, data.frame(x, y, paste0(base.url, zoom, "/", x,"/", y, ".png"), stringsAsFactors = FALSE))
	}
	}
	tiles.per.pand[[length(tiles.per.pand) + 1]] <- c(panden$pand.ID[[k]], list(l))
	}

	tiles <- unique(tiles)
	colnames(tiles) = c("x", "y", "url")
	```

	### Stap 5: Download de Zonatlas tegels en zet kleuren om in waarden
	Voordat je de Zonatlas tegels kunt downloaden en opslaan voor verdere analyse, moet je ze eerst georefereren. Hiervoor gebruik ik een aantal functies afkomstig uit de OpenStreetMap Wiki. Waar nodig heb ik de code aangepast op het tiling schema van de Zonatlas.

	```{r}
	tile2lon <- function(x, zoom) {
	return (x / 2.0 ^ zoom * 360.0 - 180.0)
	}

	tile2lat <- function(y, zoom) {
	y.OSM <- 262143 - y # LET OP: correctie t.o.v. OSM tile schema!
	n <- 2.0 ^ zoom
	lat_rad <- atan(sinh(pi * (1 - 2 * y.OSM / n)))
	lat_deg <- 180.0 * (lat_rad / pi)
	return (lat_deg)
	}

	get.bounding.box <- function(x, y, zoom){
	return (c(ymin = tile2lat(y - 1, zoom), ymax = tile2lat(y, zoom), xmin = tile2lon(x, zoom), xmax = tile2lon(x + 1, zoom)))
	}
	```

	Nu kunnen de Zonatlas tegels gedownload worden. Het zijn PNG-bestanden, die worden ingelezen en opgeslagen in het native raster formaat. Dit geeft me de mogelijkheid om kleuren in het rasterbestand om te zetten naar waarden (rood = minder geschikt, geel = geschikt, groen = zeer geschikt).

	```{r, message = FALSE, warning = FALSE, results = 'hide'}
	for (i in 1:nrow(tiles)) {
	url <- tiles$url[i]
	tmp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".png")
	download.file(url, destfile = tmp, mode = "wb")

	if (file.size(tmp) > 222) {
	img <- as.numeric(readPNG(tmp, native = TRUE))
	img[img %in% c(-14737455, -14737456)] <- 1 # rood = minder geschikt
	img[img == -14423829] <- 2 # geel = geschikt
	img[img == -14114800] <- 3 # groen = zeer geschikt
	img[!img %in% c(1, 2, 3)] <- NA
	bbox <- get.bounding.box(tiles$x[i], tiles$y[i], zoom)
	r <- raster(matrix(img, nrow = 256, ncol = 256, byrow = TRUE),
	xmn = bbox["xmin"], xmx = bbox["xmax"], ymn = bbox["ymin"], ymx = bbox["ymax"],
	crs = "+init=epsg:4326")
	writeRaster(r, paste0(tiles$x[i], "-", tiles$y[i], ".grd"), format = "raster", overwrite = TRUE)
	}
	unlink(tmp)
	}
	```

	### Stap 6: Bepaal voor ieder pand de waarde in de Zonatlas
	In de laatste stap voer ik per pand de volgende deelstappen uit:
	1. Als er meerdere Zonatlas tegels bij het pand horen, voeg de native rasterbestanden dan samen tot ��n.
	2. Selecteer alleen de rasterwaarden die binnen de begrenzing van het pand vallen.
	3. Bereken de mediaan van de rasterwaarden binnen de begrenzing van het pand.
	4. Decodeer de mediaan naar een waarde voor de geschiktheid van het pand voor het opwekken van zonne-energie

	Let op: de waarde NA geeft aan dat de geschiktheid van het pand voor het opwekken van zonne-energie niet bepaald kon worden.

	```{r, message = FALSE, warning = FALSE, results = 'hide'}
	decode.color <- function (x){
	if (is.na(x)) return (NA)
	else if (x == 1) return ("minder geschikt")
	else if (x == 2) return ("geschikt")
	else if (x == 3) return ("zeer geschikt")
	}

	for (i in 1:length(tiles.per.pand)) {
	tiles <- tiles.per.pand[[i]][[2]]
	tiles <- tiles[sapply(tiles, file.exists)]
	if (length(tiles) == 0) {
	# Voor deze panden bevat de Zonatlas geen informatie.
	panden[tiles.per.pand[[i]][[1]], c("zonenerg")] <- NA
	} else {
	if (length(tiles) == 1) {
	r <- raster(tiles[[1]])
	} else if (length(tiles) > 1) {
	r <- do.call(merge, lapply(tiles, function(x){t <- raster(x); origin(t) <- c(0,0); return (t) }))
	}
	sp <- panden[tiles.per.pand[[i]][[1]],]
	panden[tiles.per.pand[[i]][[1]], c("zonenerg")] <- decode.color(summary(mask(r, sp))[3])
	}
	}
	```

	Het eindresultaat schrijf ik tenslotte weg naar een shape-bestand met de naam panden-provincie.
	In een shape-bestand mogen attribuutnamen niet langer zijn dan 8 karakters, vandaar de lelijke namen.
	```{r, message = FALSE, warning = FALSE, results = 'hide'}
	panden <- spTransform(panden[, c("pand.ID", "zonenerg")], CRS("+init=epsg:28992"))
	sapply(c("panden-provincie.dbf", "panden-provincie.prj", "panden-provincie.shp", "panden-provincie.shx"), unlink)
	writeOGR(panden, dsn = getwd(), layer ='panden-provincie', driver = 'ESRI Shapefile')
	```

	## Resultaten
	Van de `r nrow(rijen)` objecten in de spreadsheet Gebouwen.xls uit het projecten portaal, zijn `r length(identificatienrs)` meegenomen in de analyse.
	Voor `r (nrow(rijen) - length(identificatienrs))` kon geen vlakgeometrie bepaald worden, vandaar dat ze niet meedoen.

	kable(head(panden), caption = "De eerste 6 rijen", row.names = FALSE, col.names = c("Pand ID", "Geschiktheid voor opwekken zonne-energie"), padding = 2)
	table(panden$zonenerg, useNA = "ifany")

	panden <- readOGR(dsn = getwd(), layer ='panden-provincie')
	for (i in 1:nrow(panden)) {
	panden$oppervlakte.m2[i] <- panden@polygons[[i]]@area
	}

	x <- arrange(panden@data[panden$zonenerg == 'zeer geschikt',], desc(oppervlakte.m2))
	x <- rename(x, c("pand_ID" = "BAGID"))
	data <- join(x, rijen, by = "BAGID")
	head(data[, c("BAGID", "Complex", "zonenerg", "oppervlakte.m2")], n = 5)