1

我有一个 202 行和 201 列的栅格地图 在这个地图中有一些像素值为 0 的网格 我想编写一个函数来返回所有像素值 0 网格的坐标 我该怎么做 我试图使用 if loop 和 while 循环，但它总是说 TRUE/FALSE 需要这里是我的示例代码

library(raster)
library(rgdal)
library(maptools)
library(sp)


setwd("E:/Landsat-data-NASA atm-corrected/sample_day1")
restdir2 <- ("E:/Landsat-data-NASA atm-corrected/sample_day1")
  n3 <- list.files(restdir2, pattern="*band4_clip_1.tif", full.names=TRUE)
  n4 <- list.files(restdir2, pattern="*cloud_qa_clip_1.tif", full.names=TRUE)
  n5 <- list.files(restdir2, pattern="*cloud.tif", full.names=TRUE)

create<- function(x,y)
{
 layer <- raster(n4)
 layer2 <- raster(n3)
   for(c in 1:x)
  {
    for(r in 1:y)
   {      
       nl<- layer2
       if(layer[c,r]==0)
       return layer[c,r]
   }
  }
}
create (10,10)
4

1 回答 1

1

这是两种（非常相似的）方法

library(raster)
# set up example data
r <- raster(nrow=18, ncol=36)
set.seed(0)
r[] <- round(runif(ncell(r)) * 10 - 5)

# approach 1, for a single layer
p <- rasterToPoints(r, fun=function(x){x == 0})

# approach 2, also works for multiple layers
# first remove all non zero cells
z <- subs(r, data.frame(0, 1))
p <- rasterToPoints(z)

# results
plot(r)
points(p[,1:2])

如果您有多个具有相同空间参数（范围和分辨率）的图层

# create example data 
x1 <- setValues(r, round(runif(ncell(r)) * 10 - 5))
x2 <- setValues(r, round(runif(ncell(r)) * 10 - 5))
x3 <- setValues(r, round(runif(ncell(r)) * 10 - 5))

# combine layers
s <- stack(x1, x2, x3)

z <- subs(r, data.frame(0, 1))
p <- rasterToPoints(z)
于 2015-04-21T00:52:49.773 回答