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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 |
## 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.
1 2 3 4 5 6 7 8 9 10 11 12 |
# 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) |
1 2 3 |
# 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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
# 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)
Hei!
Postnummerlistene er tilgjengelege i ei rekkje ulike format, ikkje berre HTML, men CSV, TXT, XML og KML.
Sjekk ut http://www.erikbolstad.no/geo/noreg/postnummer/last-ned/.
Ha ein fin dag!
Erik
Takk for at du no la ut listene i andre format også. Eg har oppdatert artikkelen til heller å bruka desse direkte.
Det er for øvrig ein liten feil i filformatet. CSV-fila er ikkje ei ekte CSV-fil; ho er ei TSV-fil. Skildringa «Tabseparert CSV» gjev ikkje heilt meining, då C-en i CSV nettopp står for «comma». Sjølv om nokre folk «misbrukar» namnet CSV til å tyda DSV, bør ein unngå dette. Det fører til forvirring.
Eg valde å kalle fila CSV fordi Excel og mange andre program ironisk nok ventar seg tabulatorseparerte filer framfor kommaseparerte filer når filendinga er
.csv
.No har eg endeleg lært meg å generere filer, så det er berre å kome med ynske om andre format. 🙂
Jeg kom hit, så sa det stopp.
> # Vis kor mange færre «postnummer» me no fekk
> nrow(pnr)
[1] 4665
> nrow(pnr.unik)
[1] 4331
>
>
> # 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
Error in proj4string(norge) :
error in evaluating the argument 'obj' in selecting a method for function 'proj4string': Error: object 'norge' not found
Dette tyder på at
norge
-objektet ikkje er definert hos deg. Det vert oppretta tidligare, i linjanorge=gUnaryUnion(fylke)
. Sjå om du fekk ein feilmelding når du kjørte den linja.Du har ikke nevnt at en må laste opp alle pakker du referer til, de var ikke med på standard innstallasjon.
Jeg skjønte ikke syntaksen på
fylke=readOGR(«N5000 shape», «N5000_ArealdekkeFlate»)
men fant ut etterhvert at første felt er katalogen til fila. Da gikk alt bra tilslutt.
Det jeg primært trengte var å generer en CSV-fil med UTM-xy koordinater til postnummer. Det gikk etter litt plundring.
Takk for svar!
Installering av pakkar som vert kalla er ein «standard» prosedyre i R, og derfor tok eg ikkje med informasjon om det. ☺
For berre å gjera om koordinatane til UTM-koordinatar treng du ikkje lasta inn kartet i det heile. Berre bruk denne
proj4string
-en ispTransform()
-funksjonen:CRS("+init=epsg:32633")
. Dette er koden for UTM sone 33, som ofte er brukt for å representera heile Norge. Byt ut dei to siste siffera for andre sonar. For delar av Norge er andre projeksjonar betre, frå UTM sone 31 i vest til UTM sone 36 heilt aust (i Finnmark), pluss andre projeksjonar for endå meir lokale forhold.Dette ble for vanskelig for en grønnskolling. Jeg prøvde meg på noe sånt:
> ## 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
>
>
>
> # 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=", "))
>
>
> uproj=CRS("+init=epsg:32633") # Uprojiserte data
>
> proj4string(pnr.unik)=uproj
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "proj4string pnr.proj=spTransform(pnr.unik,uproj) # Transformer til UTM sone 33
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "spTransform", for signature "data.frame", "CRS"
>
Det er fleire rare ting her. For det første bør du vel ikkje slå saman postnummer med same posisjon. For det andre har du ikkje gjort dataramma om til eit geografisk objekt (lettast med
coordinates()
-funksjonen), og då kan du ikkje gjera geografiske operasjonar på ho. For det tredje prøver du å seia at koordinatane er i UTM sone 32 (EPSG 32633), men dei er jo uprojisert (EPSG 4326), og du skal projisera dei til UTM sone 32. Og feilmeldinga du tok med ser ut til å mangla noko; ho ser iallfall ikkje rett ut.Hei,
Statenskartverk har endret/flyttet grunndata litt – jeg søkte etter N5000 og så lastet ned denne:
http://www.kartverket.no/Documents/Kart/N50-N5000%20Kartdata/33_N5000_shape.zip
og da brukte filen NO_Arealdekke_pol.shp fra zip filen.
Skal egentlig ha ut koordinater av Voronoi-tesselering for å laste inn i et annet program (Tableau) – men nå er jeg et godt stykk videre – takk for hjelpen!
For de som har tilgang til GAB (Alle adresser i Norge inkl. postnummer), så kan man lage et mer nøyaktig kart.
Først bruker man Voronoi-tesselering på alle adressene, og deretter slår man sammen alle polygonene som har samme postnummer (Dissolve). Det tar tid, men resultatet blir et helt nøyaktig postnummerkart.
Evt. kan «Digitale Postnummergrenser» (DPG) kjøpes hos Bring/Posten.