[GIS] Mask raster function in R not working to extract values in one raster based on conditions in another

extract-by-maskrraster

require: raster package

I am having trouble masking values of one raster using values in another. This is similar to this question:

https://stackoverflow.com/questions/51716315/how-to-subset-classify-raster-based-on-another-raster-grid-cells-values

but it doesn't have the moving window (focal function) as a component of the question. I have two rasters with the same dimension, extent, projection.

Raster "MAT2resampled" is a raster with a wide range of temperature indices

Raster "wintPCP" is a second raster with pixel values of 1, O, and NA

The only pixels I want from MAT2resampled are those with the value 1 in wintPCP.

I used this code:
PthreshWint2 <- mask(MAT2resampled, wintPCP)

When I plot PthreshWint2, it doesn't appear as though the MAT2resampled raster has been masked at all. I think I'm missing something simple.

I have posted the two rasters here:
https://drive.google.com/open?id=118kCu6QpoXSYmA17PpODRhHDJefHPkx1

https://drive.google.com/open?id=1N062TiHTHBqRYrclRW8O87WY-MX8taQg

Best Answer

You have a mask with 0, 1, and NA values and you want to keep only the 1s. Hence you need to convert 0 and NA to the same value.

Here's a bit of a trick. Diving m by itself results in NA where m is 0 or NA (since 0/0 returns NaN, which acts like NA in mask) and 1 where its 1:

rmask = mask(r,m/m)

not sure if that's any quicker than the other solutions.

Related Solutions

[GIS] Mosaicking rasters in R

The issue here is that mosaic and do.call are expecting a raster object in the list and not just character names of the raster that is contained in the "rasters1" vector. You are, in effect, asking to mosaic a name in a vector and not a raster object.

# Create some example data
require(raster)
    r <- raster(ncol=100, nrow=100)
      r1 <- crop(r, extent(-10, 11, -10, 11))
        r1[] <- 1:ncell(r1)
          r2 <- crop(r, extent(0, 20, 0, 20))
          r2[] <- 1:ncell(r2)
      r3 <- crop(r, extent(9, 30, 9, 30))
    r3[] <- 1:ncell(r3)

# If I create a list object of the raster names, as your are doing with list.files, 
#    do.call will fail with a character signature error 
rast.list <- list("r1","r2","r3")   
  rast.list$fun <- mean     
    rast.mosaic <- do.call(mosaic,rast.list)

# However, if I create a list contaning raster objects, the do.call function 
#   will work when mosaic is passed to it.      
rast.list <- list(r1, r2, r3)     
  rast.list$fun <- mean
    rast.mosaic <- do.call(mosaic,rast.list)
      plot(rast.mosaic)

# You could specify a for loop to create a list object, 
#   contaning raster objects
rasters1 <- list.files("F:/MOD15A2_LAI_1km/MOD15A2_LAI_2009", 
                       pattern="mod15a2.a2009001.*.005.*.img$", 
                       full.names=TRUE, recursive=TRUE)
rast.list <- list()
  for(i in 1:length(rasters1)) { rast.list[i] <- raster(rasters1[i]) }

# And then use do.call on the list of raster objects
rast.list$fun <- mean
  rast.mosaic <- do.call(mosaic,rast.list)
    plot(rast.mosaic)

[GIS] Extract the edge between two raster cells with different values

I think it is not wise to mix the clump(..) functionality from igraph with the dissolve=TRUE parameter from the rasterToPolygon routine. They both do something with to aggregate the fields together but in a different way. At least we want to do 3 things:

read or desing a raster
select raster area where the contour goes around
define how the contour is shaped (rectangular,linear or smooth) and define which areas belong together.

In the clump code raster, the type of raster and the definition of NA seems to be important to steer clumping the process. I made some test but with bad results. I followed your sketch and here is a little analysis how to get things work:

# Load packages
require('raster')
require('rgeos')

# Clean up everything
rm(list=ls())

# Set a defined random seed
set.seed(2)

# Create a float raster with values
# the interval [0,2] float
rs <- raster(nrow=10, ncol=10)

# Scale the random numbers from interval [0,1) to
# to [0,2.2) and shift the interval to [-0.1,2.1)
rs[] <- runif(ncell(rs)) * 2.2 - 0.1

# Cut off the raster values to [0,2]
# Everything smaller then 0 is zero
values(rs)[values(rs) < 0.0] <- 0.0

# Everything larger then 2 is two
values(rs)[values(rs) > 2.0] <- 2.0

Is the raster field well constructed?

# Construction is OK?
> quantile(values(rs))
       0%       25%       50%       75%      100% 
0.0000000 0.3928834 0.8733250 1.6027794 2.0000000

Here is the function prototype that selects the field in the interval [1,2).

# Function of the contour ID
inOne <- function(x) { x>=1 & x<2 }

> inOne(0) 
[1] FALSE

> inOne(1) 
[1] TRUE

> inOne(2)
[1] FALSE

The contour cannot be dissolved (clump) because of the float number nature of the raster field is distinct in the ID process (dissolve-=TRUE).

# Contour of the float desing
# x := [1,2)
ct <- rasterToPolygons(rs, 
         fun=inOne , 
         dissolve=TRUE) 
plot(rs)
plot(ct, add=TRUE)

You see the right contour groups, but polygons are not joint.

So if we have a unique ID of each cell as Integer, the dissolve process should work.

# Apply integer operation (ceiling, floor, round)
# to the float number fields 
rs.int <- ceiling(rs)
values(rs.int)
ct.int <- rasterToPolygons(rs.int, 
              fun=inOne , 
              dissolve=TRUE) 
plot(rs.int)
plot(ct.int, add=TRUE)

Conclusion: I think (do know not exactly) what behind the clump stuff works a raster based region growing routine. The process dissolve=TRUE in the rasterToPolygon (based on rgeos CRAN) seems to follow a vector approach. So I've to read the manuals of igraph and rgeos carfully.

REM: The selection of contours (float vs. int) differs, because of the nature of (ceiling, floor and round).

Best Answer

Related Solutions

[GIS] Mosaicking rasters in R

[GIS] Extract the edge between two raster cells with different values

Related Question