Karl Ove Hufthammer

Generering av områdekart basert på punktobservasjonar

27. January 2012 (oppdatert 10. april 2013)

Av og til har ein behov for å laga kart som viser korleis eit geografisk område er delt inn, basert på observasjonar i området. Denne artikkelen viset eit eksempel på korleis me lett kan gjera dette automatisk, ved av hjelp av frie verktøy.

Me skal prøva å laga eit postnummerkart, altså eit kart som viser korleis landet er delt inn i postnummersoner. Overraskande nok er ikkje dette offentlig tilgjengelig informasjon! Heldigvis har ein gjeng, med Erik Bolstad i spissen, på dugnad laga ei oversikt over midtpunktet til kvart postnummer. Denne er tilgjengelig under ein Creative Commons Navngivelse 3.0 Unported-lisens.

Me startar med å lasta ned CSV-versjonen av postnummeroversikta. Det kan òg vera greitt å bruka eit norgeskart som utgangspunkt for kartet vårt, for å hindra at me får postnummersoner som går langt til havs. Her kan me lasta ned eit gratiskart frå Statens kartverk, på litt strengare vilkår (Navngivelse-DelPåSammeVilkår 3.0 Norge). Sluttresultatet vert altså tilgjengelig under denne lisensen. Last ned shapefile-versjonen av fylkeskartet.

Til slutt treng me programpakka R, som kan brukast til det meste. For å redigera og kjøra R-skript anbefaler eg å bruka RStudio.

Så er det berre å setta i gang. Me plasserer CSV-fila med postnummer i ei mappe, og pakkar ut kartfila same plass. Litt enkel R-kode hentar inn og gjer klar dei to datakjeldene:

## Estimer postnummersoner basert på koordinatar til postnummer
 
library(rgdal)        # For lasting av kartfiler
library(rgeos)        # For geometriske operasjonar
library(plyr)         # For enkel datamanipulering
library(spatstat)     # For Voronoi-tesselering og formatkonvertering
library(spdep)        # For utrekning av polygon-naboskap
library(maptools)     # For formatkonvertering og ymse
library(RColorBrewer) # For fine fargepalettar
 
# Les inn fylkesfila
fylke=readOGR("N5000 shape", "N5000_ArealdekkeFlate")
 
# Fjern havpolygon og slå saman alle fylka til eitt stort polygon
fylke=fylke[fylke$OBJTYPE!="Havflate",]
norge=gUnaryUnion(fylke)
 
# Last inn postnummerdatabasen. Last først alt som tekst
# (for å bevara leiane nullar i postnummer), og gjer så
# lengde- og breiddegradar til tal.
pnr=read.delim("postnummer.csv", colClasses="character")
pnr=transform(pnr, LON=as.numeric(LON), LAT=as.numeric(LAT))
 
# Fjern postnummer som ikkje er i bruk
pnr.ok=pnr[!grepl("ikkje i bruk", pnr$POSTSTAD),]
 
# Nokre posisjonar er dei same for ulike postnummer. Dette
# fører til problem, så slå dei derfor saman til eitt nummer.
pnr.unik=ddply(pnr.ok, ~LON+LAT, summarise, POSTNR=paste(POSTNR, collapse=", "))
 
# Vis kor mange færre «postnummer» me no fekk
nrow(pnr)
nrow(pnr.unik)
 
 
# Norgeskartet er i UTM-projeksjon sone 33, mens postnummera 
# er uprojiserte (dvs. lengde- og breiddegradar). Gjer derfor 
# koordinatane til postnummera om til UTM sone 33-koordinatar.
proj=CRS(proj4string(norge)) # UTM sone 33
uproj=CRS("+init=epsg:4326") # Uprojiserte data
 
coordinates(pnr.unik)=~LON+LAT       # Gjer om til geografiske objekt
proj4string(pnr.unik)=uproj          # Legg til info om projeksjon
pnr.proj=spTransform(pnr.unik, proj) # Transformer til UTM sone 33
 
# Sidan norgeskartet vårt berre dekker fastlandet, kutt
# ut postnummer på Svalbard og deromkring, samt postnummer
# som ikkje er i bruk
pnr.ok=pnr.proj[!is.na(over(pnr.proj, norge)),]
 
# Endå færre postnummer no ...
nrow(pnr.ok)
 
# Vis postnummera på norgeskartet
plot(norge, border="darkgrey")
points(pnr.ok, col="red", pch=".")

For å «gjetta» oss fram til korleis kvar postnummersone ser ut, lagar me ei Voronoi-tesselering. Wikipedia har ei veldig fin forklaring på kva dette er, men kort sagt går det ut på at me for kvar postnummerposisjon generer ei postnummersone som består av alle geografiske punkta som er nærare denne posisjonen enn dei er nokon av dei andre postnummerposisjonane. Ut av dette får me ei postnummer-tesselering (flisinndeling) som dekker heile landet.

# For å laga Voronoi-tesseleringa går me ein liten omveg
# og brukar «spatstat»-pakken, sidan me lett kan gjera
# resultatet om til eit geografisk objekt.
bboks.norge=bbox(norge)
bboks.owin=owin(bboks.norge[1,], bboks.norge[2,])
pnr.ppp=as.ppp(coordinates(pnr.ok), W=bboks.owin)
fliser=dirichlet(pnr.ppp)
fliser.poly=as(fliser, "SpatialPolygons")
proj4string(fliser.poly)=proj
 
# Resultatet er ei fin tesselering
plot(fliser.poly)

# Fjerner så alt som ligg utanfor Norges landegrenser
fliser.poly.norge=gIntersection(fliser.poly, norge, byid=TRUE)
plot(fliser.poly.norge)

Dette ser jo greitt ut, men det kan vera fint med litt fargar òg. Som kjent finst det eit matematisk teorem som seier at alle kart kan fargeleggast med maks fire fargar slik at ingen naboland har same farge. Nokre kart kan vera fryktelig vanskelige (men mulig!) å farga med fire fargar, så me nøyer oss med ein enklare og kjapp algoritme, som likevel brukar få fargar. Det finst postnummersoner med både 12 og 13 nabosoner, men denne algoritmen klarar å fargelegga kartet med berre seks (nesten fem!) fargar.

# Lag fargetal/-indeksar for fargelegging av kart
lagKartFargar=function(x) {
  nb=poly2nb(x)          # Lag naboliste
  n=length(x)            # Talet på polygon
  cols=numeric(n)        # Lag klar fargevektor
  cols[1]=1              # La første polygon ha farge 1
  cols1n=1:n             # Tilgjengelige fargeindeksar
  for(i in 2:n)          # For kvart polygon, bruk «lågaste» tilgjengelige farge
    cols[i]=which.min(cols1n %in% cols[nb[[i]]])
  cols
}
 
# For å bruka minst mulig fargar, kan det hjelpa
# å sortera polygona, slik at me først prøver
# å fargelegga polygonet med fleste naboar, så
# det mest nest mest naboar, osv.
naboar = poly2nb(fliser.poly.norge, queen=FALSE)
ord = order(card(naboar), decreasing=TRUE)
fliser.poly.norge = fliser.poly.norge[ord,]
 
cols = lagKartFargar(fliser.poly.norge) # Rekn ut fargeindeksar
ncols = max(cols)                       # Talet på fargar som vart brukt
ncols                                   # Sjå på talet på fargar
 
# Grei fargepalett (prøv «display.brewer.all()» for alternativ)
pal = brewer.pal(ncols, "Set1")
 
# Teikn opp postnummerkartet
plot(fliser.poly.norge, col=pal[cols], border="transparent")

Så enkelt var det å laga eit postnummerkart! Og når me først har dataa, kan me eksperimentera litt. Kva med eit kart over første to siffer i postnummera?

Samanlikn gjerne dette med Posten Norge sitt kart (legg for øvrig merke til den ikkje heilt optimale fargelegginga der). Og hugs at vårt kart er basert på automatisk inndeling i postnummersoner, utan tilgang til dei faktiske grensene.

Me kan òg laga eit kart over berre første siffer i postnummera.

Det ser litt ut som Norge etter mulige framtidige fylkessamanslåingar! ☺

Til slutt: Korleis kunne postnummerkartet gjerast betre? Jo, viss me hadde fleire posisjonar frå kvar postnummersone, kunne me fått mykje meir nøyaktige sonegrenser. Hugs at i vårt kart er kvar sone basert på berre éin posisjon (pluss indirekte informasjon frå nabopostnummer)!

Oppgåva til lesaren som vil øva seg litt på handtering av slike geografiske data: Last ned heile programmet ovanfor. Test korleis fleire observasjonar innanfor kvar sone påverkar kor nøyaktig sluttkartet vert. Som utgangspunkt kan du prøva å generera eit fylkeskart berre basert på tilfeldige posisjonar i ulike fylke. Bruk funksjonen spsample() til å plukka nokre posisisjonar frå kartobjektet fylke. (Dette er i prinsippet omtrent det same som å ringa til ein haug tilfeldige personar i Norge og be dei oppgje geografisk posisisjon og fylke.) Kjør så ei tesselering. Kvar observasjon vert òg her til éi sone. Slå saman soner i same fylke ved å bruka unionSpatialPolygons()-funksjonen. Samanlikn til slutt grafisk dei automatisk genererte fylkesgrensene med «fasiten» i fylke-objektet.

Kjelder:
Postnummerliste: www.erikbolstad.no/postnummer/
Kartgrunnlag: Statens kartverk (cc-by-sa-3.0)

Emne: Geodata, Matematikk og statistikk, Programvare

[Abonner på kommentarar til artikkelen]

Kommentarar

Legg til kommentar





Du kan bruka dei vanligaste elementa og attributta i HTML. Avsnitt lagar du med vanlige linjeskift. Eg kan komma til å gjera typografiske og ortografiske endringar i innlegg (men vil aldri endra sjølve innhaldet), samt fjerna upassande innlegg.

Skriven av Karl Ove Hufthammer og driven med WordPress. Du kan abonnera på innleggs-RSS eller kommentar-RSS.