alternative bivariate chloropleth map legend done in R with ggplot

alternativeBivariateLegend.R

# alternative bivariate chloropleth map legend
# by Florian Hoedt
# gis.hoedt@gmail.com

require(maptools)
require(ape) #needed?

# shp with data
dataPath <- "C:/temp/Modul_09_VisKarto/A/A_1/both.shp"
shpData <- readShapePoly(dataPath)

### Fertility
f <- shpData$FertRate
f.noNull <- subset(f,!is.na(f))
f.max <- max(f.noNull)

hist(f,
     main = "Geburtenrate bei Jugendlichen",
     breaks=8,
     xlim=c(0, f.max * 1.2),
     #ylim=c(0,100),
     xlab = "Im Alter von 15-19 Jahren [Geburt je 1000 Frauen]",
     sub = "blaue Linien entsprechen Klassengrenzen",
     col="gray",
     border="white")
rug(f,ticksize = -0.05,col="red")
f.groups = c(88,121)
for (i in f.groups)
{
  abline(v=i,lty=2,col="blue")
}

### AGRICULTURE
a <- shpData$AgrRate
a.noNull <- subset(a,!is.na(a))
a.max <- max(a.noNull)

hist(a,
     main = "Landwirtschaftlicher Flächenanteil",
     breaks=8,
     xlim=c(0, a.max * 1.2),
     #ylim=c(0,100),
     xlab = "Im Alter von 15-19 Jahren [Geburt je 1000 Frauen]",
     sub = "blaue Linien entsprechen Klassengrenzen",
     col="gray",
     border="white")
rug(a,ticksize = -0.05,col="red")
a.groups = c(42,65)
for (i in a.groups)
{
  abline(v=i,lty=2,col="blue")
}

#Farben
b.A3.col = rgb(108,131,181, maxColorValue=255)
b.B3.col = rgb(88,122,149, maxColorValue=255)
b.C3.col = rgb(41,91,92, maxColorValue=255)

b.A2.col = rgb(180,192,218, maxColorValue=255)
b.B2.col = rgb(144,178,179, maxColorValue=255)
b.C2.col = rgb(91,145,121, maxColorValue=255)

b.A1.col = rgb(232,232,232, maxColorValue=255)
b.B1.col = rgb(184,214,190, maxColorValue=255)
b.C1.col = rgb(115,174,130, maxColorValue=255)

#bins
b.Al <- 0
b.Ah <- 41
b.Bl <- 41
b.Bh <- 64
b.Cl <- 64
b.Ch <- 80

b.1l <- 0
b.1h <- 88
b.2l <- 88
b.2h <- 120
b.3l <- 120
b.3h <- 210

# Leerplot (type="n") im Koordinatenbereich der st 
plot(x=a, y=f, type="n",
     ylab="Geburtenrate bei Jugendlichen [Geburt je 1000 Frauen]", xlab="Landwirtschaftlicher Flächenanteil [% Gesamtfläche]")
# Einfügen der Teilgebietsgrenzen
#B.A1
rect( ytop    = b.1h, 
      xleft   = b.Ah,
      xright  = b.Al,
      ybottom = b.1l,
      border = "dark gray", col = b.A1.col)

#B.A2
rect( ytop    = b.2h, 
      xleft   = b.Ah,
      xright  = b.Al,
      ybottom = b.2l,
      border = "dark gray", col = b.A2.col)

#B.A3
rect( ytop    = b.3h, 
      xleft   = b.Ah,
      xright  = b.Al,
      ybottom = b.3l,
      border = "dark gray", col = b.A3.col)

#B.B1
rect( ytop    = b.1h, 
      xleft   = b.Bh,
      xright  = b.Bl,
      ybottom = b.1l,
      border = "dark gray", col = b.B1.col)

#B.B2
rect( ytop    = b.2h, 
      xleft   = b.Bh,
      xright  = b.Bl,
      ybottom = b.2l,
      border = "dark gray", col = b.B2.col)

#B.B3
rect( ytop    = b.3h, 
      xleft   = b.Bh,
      xright  = b.Bl,
      ybottom = b.3l,
      border = "dark gray", col = b.B3.col)

#B.C1
rect( ytop    = b.1h, 
      xleft   = b.Ch,
      xright  = b.Cl,
      ybottom = b.1l,
      border = "dark gray", col = b.C1.col)

#B.C2
rect( ytop    = b.2h, 
      xleft   = b.Ch,
      xright  = b.Cl,
      ybottom = b.2l,
      border = "dark gray", col = b.C2.col)

#B.C3
rect( ytop    = b.3h, 
      xleft   = b.Ch,
      xright  = b.Cl,
      ybottom = b.3l,
      border = "dark gray", col = b.C3.col)


# für cor() dürfen keine NA Werte enthalten sein
# z.B. Definitionsabfrage in ArcGIS nutzen und neuen bereinigten
# Datensatz exportieren
dataPath <- "C:/temp/Modul_09_VisKarto/A/A_1/bothNoNull.shp"
shpDataNoNull <- readShapePoly(dataPath)
corrFA <- cor(shpDataNoNull$AgrRate,shpDataNoNull$FertRate)

# Einfügen der Baumpunkte
points(x=a,y=f,pch = ".", cex=5,
       col="black")
# Einfügen Titel und Untertitel
title("Die echte Aufteilung der Klassen",
      sub = paste0("Korrelationskoeffizient Landwirtschaftliche Flächen zu Geburtenrate: ",round(corrFA,digits = 4)),
      cex.sub = 0.75, font.sub = 3, col.sub = "dark gray")