2013-11-28 40 views
0

我被分配,我想在R.尝试shape文件(SpatialPoints DF)被导入由几个属性为特定的任务,但最重要的是,商业捕捞权重具体点坐标(纬度/经度)。汇总样本由方格

我要求的脚本:

1)创建了一个网格(其中尺寸和单元可以修改) 2)与进口文件相交为了总结样品(平均,SD,范围等)由网格平方。

我可以通过ArcGIS这样做,但很感兴趣,便于修改网格大小和具有通过R.可重复使用的算法的下面是所使用的数据的一个简短的例子。

会有人知道如何做到这一点?

ENT_LATITU ENT_LONGIT CSK 
    415300  654400 195.954 
    430100  622200 21.228 
    442300  631000 232.423 
    424700  642300 77.837 
    442800  630600 154.586 
    424600  642900 9.253 
+0

这将是更好,如果你提供创建对象的shp-码对象和类似的范围规格您将与之合作的对象。 –

回答

2

我建议使用raster & sp软件包来汇总网格单元格中的点。下面的代码应该让你开始。这允许您设置行数和列数,如果你想要设置单元格的大小,而不是它不会是很难对此进行修改。

library(sp) 
library(raster) 

#recreate a sample of your data 
dF <- data.frame(ENT_LATITU=c(415300,430100,442300,424700),ENT_LONGIT=c(654400,622200,631000,642300), CSK=c(195,21,232,77)) 

nameLon <- "ENT_LATITU" 
nameLat <- "ENT_LONGIT" 

#put points into a 'SpatialPointsDataFrame' 'sp' object 
coords <- cbind(x=dF[[nameLon]],y=dF[[nameLat]]) 
sPDF <- SpatialPointsDataFrame(coords,data=dF) 

#set number of rows & columns in the grid 
nrows <- 3 
ncols <- 3 

#setting extents from the data 
xmn <- min(dF[[nameLon]]) 
ymn <- min(dF[[nameLat]]) 
xmx <- max(dF[[nameLon]]) 
ymx <- max(dF[[nameLat]]) 

#create a grid 
blankRaster <- raster(nrows=nrows, ncols=ncols, xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx) 
#adding data into raster to avoid 'no data' error 
blankRaster[] <- 1:ncell(blankRaster) 

#calc mean (or other function) of points per cell 
rasterMeanPoints <- rasterize(x=sPDF, y=blankRaster, field='CSK', fun=mean) 

#plot to get an idea whether it's doing the right thing 
plot(rasterMeanPoints) 
text(sPDF,sPDF$CSK) 
+0

我要感谢你@Andy作出的贡献。回过头来开发这个scirpt,结果我用了很多你的代码来让我通过我自己的。很多apprecaited。 – RyMar

0

如果你有一个名为“逸”到要应用功能mean通过他们的中点分裂这两套坐标形成范围的数据对象:

# dput(dat) 
dat <- 
structure(list(ENT_LATITU = c(415300L, 430100L, 442300L, 424700L, 
442800L, 424600L), ENT_LONGIT = c(654400L, 622200L, 631000L, 
642300L, 630600L, 642900L), CSK = c(195.954, 21.228, 232.423, 
77.837, 154.586, 9.253)), .Names = c("ENT_LATITU", "ENT_LONGIT", 
"CSK"), class = "data.frame", row.names = c(NA, -6L)) 

with(dat, tapply(CSK, list(lat.cut=cut(ENT_LATITU, 2), 
          lon.cut=cut(ENT_LONGIT ,2)), 
         mean)) 
#-------------------------------- 
        lon.cut 
lat.cut    (6.22e+05,6.38e+05] (6.38e+05,6.54e+05] 
    (4.15e+05,4.29e+05]     NA    94.348 
    (4.29e+05,4.43e+05]    136.079     NA 

这给了你一个表格对象(从矩阵类继承)。

1

我认为你应该使用ggplot2geom_raster()。这是一个使用一些合成数据的例子。我首先创建了一个30x30的网格,然后展示如何将其修剪为任何x/y聚合。

require(ggplot2) 
require(plyr) 

## CREATE REASONABLE SIZE GRID 30x30 
dfe<-expand.grid(ENT_LATITU=seq(415000,418000,100), 
      ENT_LONGIT=seq(630000,633000,100), 
      CSK=0) 
## FILL WITH RANDOM DATA 
dfe$CSK=round(rnorm(nrow(dfe),200,50),0) 

####################################################### 
##### VALUES TO CHANGE IN THIS BLOCK    ##### 
####################################################### 
## TRIM ORIGINAL DATASET 
lat.max<-Inf  # change items to trim data 
lat.min<-0  
long.max<-Inf  
long.min<-631000  
dfe.trim<-dfe[findInterval(dfe$ENT_LATITU,c(lat.min,lat.max))*findInterval(dfe$ENT_LONGIT,c(long.min,long.max))==1,] 
## SUMMARIZE TO NEW X/Y GRID 
xblocks<-6 
yblocks<-8 

## GRAPH COLOR AND TEXT CONTROLS 
showText<-TRUE 
txtSize<-3 
heatmap.low<-"lightgreen" 
heatmap.high<-"orangered" 
####################################################### 
#####            ##### 
####################################################### 

## BASIC PLOT (ALL DATA POINTS) 
ggplot(dfe) + 
    geom_raster(aes(ENT_LATITU,ENT_LONGIT,fill=CSK)) + theme_bw() + 
    scale_fill_gradient(low=heatmap.low, high=heatmap.high) + 
    geom_text(aes(ENT_LATITU,ENT_LONGIT,label=CSK,fontface="bold"), 
      color="black", 
      size=2.5) 

基本情节:

enter image description here

然后聚集的情节:

## CALL ddply to roll-up the data and calculate summary means, SDs,ec 
dfe.plot<-ddply(dfe.trim, 
     .(lat=cut(dfe.trim$ENT_LATITU,xblocks), 
     long=cut(dfe.trim$ENT_LONGIT,yblocks)), 
     summarize, 
     mean=mean(CSK), 
     sd=sd(CSK), 
     sum=sum(CSK), 
     range=paste(min(CSK),max(CSK),sep="-")) 

## BUILD THE SUMMARY CHART 
g<-ggplot(dfe.plot) + 
    geom_raster(aes(lat,long,fill=sum),alpha=0.75) + 
    scale_fill_gradient(low=heatmap.low, high=heatmap.high) + 
    theme_bw() + theme(axis.text.x=element_text(angle=-90)) + 
    ggtitle(paste(xblocks, 
       " X ", 
       yblocks, 
       " grid of Catch Data\nbetween (", 
       min(dfe.trim$ENT_LATITU), 
       " : ", 
       min(dfe.trim$ENT_LONGIT), 
       ") and (", 
       max(dfe.trim$ENT_LATITU), 
       " : ", 
       max(dfe.trim$ENT_LONGIT), 
       ")\n\n", 
       sep="")) 

## ADD THE LABELS IF NEEDED 
if(showText)g<-g+geom_text(aes(lat,long,label=paste("SUM=",round(sum,0), 
              "\nMEAN=",round(mean,1), 
              "\nSD=",round(sd,1), 
              "\nRNG=",range,sep=""), 
            fontface=c("italic")), 
            color="black",size=txtSize) 

## FUDGE THE LABELS TO MAKE MORE READABLE 
## REPLACE "," with newline and "]" with ")" 
g$data[,1:2]<-gsub("[,]",replacement=" to\n",x=as.matrix(g$data[,1:2])) 
g$data[,1:2]<-gsub("]",replacement=")",x=as.matrix(g$data[,1:2])) 

## PLOT THE CHART 
g + labs(x="\nLatitude", y="Longitude\n", fill="Sum\nBlock\n") 

## SHOW HEADER OF data.plot 
head(dfe.plot) 

enter image description here

+0

谢谢大家的帮助。特洛伊,你的方法就是我正在寻找的东西,情节是一个额外的好处。非常感激。 – RyMar

+0

欢迎你 - 我一直在做一些非常相似,因此很容易端口。 – Troy