r spatial mapping

Question

I am trying to make some nice spatial maps in R. I am trying to understand how to upload the data in case you would like to have a look at them but I have not figured that out (I am sorry for that but being a new user means looking for all these things).

What is my situation. I have a shapefile of the whole USA, I only need some states, and I can select my grid when I plot it (as you see from the code in the plot section).

I also have some yield data (points) which have latitude, longitude and Yield. I have 4 different Yield data which are called ("All Stations", "0.5", "1.0", and "2.0").

I am trying to plot these 4 yields data on the spatial map to have 4 different spatial maps. Which is done.

I have done this by reading on stackoverflow here and there I used bits and pieces to do that, although I have never done it before I was surprised how fast I could advance (Thank you people on StackOverflow!!!).

Can someone help me to understand if my code is correct?

Also, how can I make the legends of the 4 maps a regularr scale? E.g. from 4000 to 9000 with 500 intervals for each of the 4 maps. What I have done is to create a separate text file ("Yield for Legend.txt") which I use to generate the colour scales in the maps and the legends. Is that correct?

Again, your critics are mostly welcome!

Thank you, David

    rm(list=ls())
    setwd("C:\\Users\\.....\\Shape File")
    library(spatstat)
    library(rgdal)
    library(shapefiles)
    library(maptools)
    library(RColorBrewer)
    library(classInt)


    # read in shapefiles

    counties.rg <- readOGR("C:\\Users\\......\\Shape File", "tl_2011_us_county")


    Yields <- read.table("Yield.txt", skip=1, header = F)

    Yield.g <- as.ppp(Yields, owin( c(-89, -76), c(25, 37)))

    ## Reading Data for Legend and colouring breaks

     Y.LE <- read.table("Yield for Legend.txt", header=F)
     Y.L.I <- classIntervals(Y.LE$V1, n=9, style = "quantile")
     Y.L.I <- Y.L.I$brks

    #select color palette and the number colors (levels of income) to represent on the map
    #colors <- brewer.pal(9, "RdYlGn")
    colors <- brewer.pal(9, "Greys")

    ################################################
    ### Generating MAPS ############################
    ################################################

    #set breaks for the 9 colors
      #par(mfrow=c(2,2))

     pdf("13 August Spatial Maps.pdf") 
    # All Points 

    brks.all <-classIntervals(Yields$V3, n=9, style = "quantile")
    brks.all <- brks.all$brks

    plot(counties.rg, axes=TRUE, border="grey", xlim = c(-82, -80),
        ylim = c(24, 37))
    points(Yield.g, cex= 1.1, bg=colors[findInterval(Yields$V3, Y.L.I,all.inside=TRUE)], pch=21)

    #add a title
    title(paste ("Rainfed Yield (kg/ha)All Stations"))

    #add a legend
    legend("bottomright", legend=leglabs(round(Y.L.I)), fill=colors, bty="n", cex=0.7 ) #,x.intersp = .5, y.intersp = .5)



    # 0.5 Grid  
    brks.05 <-classIntervals(Yields$V4, n=9, style = "quantile")
    brks.05 <- brks.05$brks

    plot(counties.rg, axes=TRUE, border="grey", xlim = c(-82, -80),
        ylim = c(24, 37))
    points(Yield.g, cex= 1.1, bg=colors[findInterval(Yields$V4, Y.L.I,all.inside=TRUE)], pch=21)
    #abline(v=GF$V1, col="grey40")
    #abline(h=GF$V2, col="grey10", lty="dotted")
    #backup
    #points(Yield.g, cex= Yields$V4/9000, col=colors[findInterval(Yields$V4, brks.05,all.inside=TRUE)], pch=19)

    #add a title
    title(paste ("Rainfed Yield (kg/ha)0.5"))

    #add a legend
    legend("bottomright", legend=leglabs(round(Y.L.I)), fill=colors, bty="n", cex=0.7 ) #,x.intersp = .5, y.intersp = .5)


    # 1.0 Grid
    brks.1 <-classIntervals(Yields$V5, n=9, style = "quantile")
    brks.1 <- brks.1$brks

    plot(counties.rg, axes=TRUE, border="grey", xlim = c(-82, -80),
        ylim = c(24, 37))
    points(Yield.g, cex= 1.1, bg=colors[findInterval(Yields$V5, Y.L.I,all.inside=TRUE)], pch=21)
    #abline(v=GO$V1, col="grey40")
    #abline(h=GO$V2, col="grey10", lty="dotted")
    #add a title
    title(paste ("Rainfed Yield (kg/ha)1.0"))

    #add a legend
    legend("bottomright", legend=leglabs(round(Y.L.I)), fill=colors, bty="n", cex=0.7 ) #,x.intersp = .5, y.intersp = .5)


    # 2.0 Grid

    brks.2 <-classIntervals(Yields$V6, n=9, style = "quantile")
    brks.2 <- brks.2$brks

    plot(counties.rg, axes=TRUE, border="grey", xlim = c(-82, -80),
        ylim = c(24, 37))
    points(Yield.g, cex= 1.1, bg=colors[findInterval(Yields$V6, Y.L.I,all.inside=TRUE)], pch=21)
    #abline(v=GG$V1, col="grey40")
    #abline(h=GG$V2, col="grey10", lty="dotted") 
    #add a title
    title(paste ("Rainfed Yield (kg/ha)2.0"))

    #add a legend
    legend("bottomright", legend=leglabs(round(Y.L.I)), fill=colors, bty="n", cex=0.7 ) #,x.intersp = .5, y.intersp = .5)
dev.off()

Answer 1

In response to your specific query about setting the breaks and having your legend match, there is a fairly straightforward fix.

You use the classInterval function with style="quantile" to define your breaks. If you want the maps to show "from 4000 to 9000 with 500 intervals for each of the 4 maps" why not use style="fixed"

brks.all <-classIntervals(Yields$V3, n=10, style = "fixed",
  fixedBreaks=seq(from=4000, to=9000, by=500)
brks.all <- brks.all$brks 

Note that 4k to 9k by 500's creates 10 intervals by my count and 9 color gradations doesn't often make for a pretty map.

Alternatively, the dataPrecision variable within classInt may also help you get labels and breaks that are closer to what you want, but still based on quantiles (if not uniform across maps)