2017-05-24 84 views
0

我想用包Spatstat創建協變量,如果點在多邊形之外,我有一個shp文件作爲多邊形以便爲協變量賦值。我的問題是在座標投影上存在問題,因爲當我完成這個過程時,我有一個奇怪的數字,請有人幫助我。我的代碼:R:在投影座標系上創建一個協變量

#Creating covariate Z: 

W<- owin(xrange=c(767768.3,768773.6),yrange=c(517335.8,518044.8)) 
xp<-c(768773.6) #cambia! 
yp<-c(518044.8) 
X<-ppp(x=xp,y=yp,W) 
par(mfrow=c(2,2)) 
Z<-density(X,500)/(7e-06) 

#reading and projecting the shp: 

Prodes1<-readOGR(dsn="forestsseparate.shp", layer="forestsseparate", dropNULLGeometries=TRUE) 
projection(Prodes1)<-CRS("+proj=utm +zone=18 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs") 

#Extracting the shp polygon coordinates: 

xcoords<-rev([email protected][[1]]@Polygons[[1]]@coords[,1]) 
ycoords<-rev([email protected][[1]]@Polygons[[1]]@coords[,2]) 
pruebalu<-owin(poly=list(x=xcoords,y=ycoords)) 

mx<-(Z$xrange[2]-Z$xrange[1])/Z$xstep 
my<-(Z$yrange[2]-Z$yrange[1])/Z$ystep 

xcoorde<-matrix(NA,nrow=1,ncol=mx) 
ycoorde<-matrix(NA,nrow=1,ncol=my) 

#Finding each pixel center 

for(i in 1:mx){ 
    xcoorde[i]<-(Z$xrange[1]+i*Z$xstep) -Z$xstep/2 
} 

for(i in 1:my){ 
    ycoorde[i]<-(Z$yrange[1]+i*Z$ystep) -Z$ystep/2 
} 

#Changing the pixel value if is outside the polygon: 

for(i in 1:length(Z$xcol)){ 
    for(j in 1: length(Z$yrow)){ 

     if(inside.owin(xcoorde[i],ycoorde[j], pruebalu)==FALSE){ 
      whichPixel<-nearest.raster.point(xcoorde[i], ycoorde[j], Z) 
      Z$v[whichPixel$col,whichPixel$row]<-0.4 
     } 
    } 
} 

這裏我的結果:

enter image description hereenter image description hereenter image description here

回答

0

對不起,我發與Z $ V族元素的順序搞錯,這裏是修正:

for(i in 1:length(Z$xcol)){ 
    for(j in 1: length(Z$yrow)){ 

     if(inside.owin(xcoorde[i],ycoorde[j], pruebalu)==FALSE){ 
      whichPixel<-nearest.raster.point(xcoorde[i], ycoorde[j], Z) 
      Z$v[whichPixel$row,whichPixel$col]<-0.4 
     } 
    } 
} 

現在完美!