1

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()
4

1 回答 1

2

针对您关于设置休息时间和让您的图例匹配的具体查询，有一个相当简单的解决方法。

您可以使用带有 style="quantile" 的 classInterval 函数来定义您的休息时间。如果您希望地图显示“从 4000 到 9000，4 个地图中的每一个都有 500 个间隔”，为什么不使用 style="fixed"

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

请注意，按我的计数，4k 到 9k 到 500 会创建 10 个间隔，而 9 个颜色渐变通常不会形成漂亮的地图。

或者，classInt 中的 dataPrecision 变量也可以帮助您获得更接近您想要的标签和中断，但仍基于分位数（如果跨地图不统一）

于 2012-08-13T19:47:45.593 回答