2016-09-21 1 views
1

Jusqu'à présent, j'utilisé la déclaration suivante dans la calculatrice raster d'ArcGIS:commande conditionnelle du transfert Raster Calculator à R

Con(("Land_use.rst" == -20), "Export.rst") 

Ce calcule une nouvelle Raster qui ne contient que les données d'exportation où Land_use est égal à -20. C'est exactement ce que je veux. Mais je veux automatiser ceci dans R, comme je dois le faire beaucoup de fois.

Jusqu'à présent, j'obtenu quelque chose comme ceci:

for (catch_dir in Dir_List) { 

r1 <- raster(paste0(catch_dir, '/Export.rst')) 

r2 <- raster(paste0(catch_dir,'/LAND_use.rst')) 

### statement that output export_streams.rst contains the data of export.rst where LAND_use.rst equals -20. 
    x <- overlay(r2, r1, fun=function(x,y){ y[x!=-20] <- 0; y}) 

writeRaster(x, filename = paste0(catch_dir, "/export_streams.rst"), format="IDRISI", NAflag = 0, datatype = "FLT4S", 
       overwrite=TRUE) 
} 

Ce code fait ce qui suit: Il produit une trame qui contient les valeurs de r1 où r2 = -20. Donc c'est bien, mais il y a aussi des restes à la frontière du raster qui ne sont pas égaux à -20 en r2. L'étendue des deux rasters est la même. Donc, la superposition ne fonctionne probablement pas pour ma tâche. D'autres idées que de superposer? Peut-être une commande if(r2 == -20){ }?

Répondre

0

Je n'ai pas vos données, donc je ne peux pas reproduire votre erreur. Cependant, voici comment je le ferais avec des données factices:

#create dummy matrix 
a <- matrix(c(10,10,10,10,10,10,10,10,10),nrow=3) 
b <- matrix(c(-20,-20,-20,3,3,3,4,4,4),nrow=3) 

#create raster from matrix 
r1 <- raster(a) 
r2 <- raster(b) 

#Create raster with same extent and nrows/ncols as original rasters 
x <- raster(extent(r1), nrows=nrow(r1), ncols=ncol(r1)) 

#produce raster x which contains the values of r1 where r2 = -20 
x[r2[]==-20] <- r1[r2[]==-20]