3

I'm running into an unusual issue with the cover function in R. We are trying to fill cloudy pixels with values from another layer. I can make it work just fine with stacks as follows -

library(raster)
r1 <- raster(ncols=36, nrows=18)
r1[] <- 1:ncell(r1)
r1b <- r1a <- r1
r1_stack <- stack(r1, r1a, r1b)

r2 <- setValues(r1, runif(ncell(r1)))
r2b <- r2a <- r2
r_stack <- stack(r2, r2a, r2b)

r_stack[r_stack < 0.5] <- NA

r3 <- cover(r_stack, r1_stack)

But then i try to do the same thing with a raster stack and i get the error:

Error in as.character(x) : 
 cannot coerce type 'closure' to vector of type 'character'

The code:

# get all tifs
LS5_032_032_2008_09_21 <- list.files("LT050340302008090301T1-SC20170526100900/", 
                                     pattern = glob2rx("*band*.tif$"), full.names = T) 


# stack bands
cloudy_scene <- stack(LS5_032_032_2008_09_21)
# import cloud mask
cloud_mask <- raster('LT050340302008090301T1-SC20170526100900/LT05_L1TP_034030_20080903_20160905_01_T1_sr_cloud_qa.tif')

# mask data 
masked_data <- mask(cloudy_scene, mask = cloud_mask, maskvalue=0, inverse=TRUE)

####### get cloud free data
# get files
LS5_2008_09_19 <- list.files("LT050340302008091901T1-SC20170526101124/", 
                             pattern = glob2rx("*band*.tif$"), full.names = T) 

# subset and stack cloud free bands
cloud_free_data <- stack(LS5_2008_09_19)

# use cover function to assign NA pixels to corresponding pixels in other scene
cover <- cover(masked_data, cloud_free_data)

TRACEBACK() output: 

9: toupper(format)
8: .defaultExtension(format)
7: .getExtension(filename, filetype)
6: .local(x, filename, ...)
5: writeStart(outRaster, filename = filename, format = format, datatype = datatype, 
       overwrite = overwrite)
4: writeStart(outRaster, filename = filename, format = format, datatype = datatype, 
       overwrite = overwrite)
3: .local(x, y, ...)
2: cover(masked_data, cloud_free_data)
1: cover(masked_data, cloud_free_data)

UPDATE: I tried to resample the data - still doesn't work

cloud_free_resam <- resample(cloud_free_data, masked_data)
cover <- cover(masked_data, cloud_free_resam)

ERROR:

Error in as.character(x) : 
  cannot coerce type 'closure' to vector of type 'character'

I also tried to crop both layers - same error 

# find intersection boundary
crop_extent <- intersect(extent(cloud_free_data), extent(masked_data))
cloud_free_data <- crop(cloud_free_data, crop_extent)
masked_data <- crop(masked_data, crop_extent)

# use cover function to assign NA pixels to corresponding pixels in other scene
cover <- cover(masked_data, cloud_free_data)

GET THE DATA: (WARNING: 317mb download - unpacks to ~1gb) https://ndownloader.figshare.com/files/8561230

Any ideas what might be causing this error with this particular dataset? I'm sure we're missing something quite basic but...what? Thank you in advance.

Leah

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
Leah Wasser
  • 717
  • 3
  • 8
  • 22
  • Works for me if I crop everything to a small extent defined by `e = extent(450000,470000,4750000,4780000)`, so eg: `cloud_free_data <- crop(stack(LS5_2008_09_19[5:10]),e)`, My PC grinds badly on the full data, I've not let it run to completion/error. – Spacedman Jun 01 '17 at 16:08
  • What does `traceback()` say after your error? Can you edit that into the question? – Spacedman Jun 01 '17 at 16:13
  • @Spacedman i have a feeling that is has something to do with overlapping pixels in areas where there are no data but i'm not certain. i'll try the crop and the traceback – Leah Wasser Jun 01 '17 at 16:30
  • @Spacedman done! plus i also uploaded a smaller dataset - it's about 1GB unpacked and 317MB zipped. – Leah Wasser Jun 01 '17 at 16:32
  • it works for me as well, when both datasets are cropped. so it could be an extent overlap issue or pixel overlap issue. – Leah Wasser Jun 01 '17 at 16:35

2 Answers2

1

I think it is to do with when the object is converted to a rasterBrick and the raster written to a temporary file. i.e., masked_data <- mask(cloudy_scene, mask = cloud_mask)

Using Spacedman's crop creates a 'in memory' rasterBrick objects that do not rely on accessing the files; in that case the example works with no error.

But using the full extent (or cropping at end of process), raster writes and accesses (temporary) files and the error occurs.

Maybe as a temporary fix is to explicitly split up the images into memory sized chunks, mask/cover and then restitch with mosaic.

The extents of the two objects are also slightly different, so that should be fixed. Also generally a good idea to avoid issues by explicitly setting band names, min-max values and values < 0 to NA.

Edward P. Morris
  • 71
  • 2
  • 7
  • Thank you! the extent issue is a problem. I tried cropping both to the same extent but that didn't fix the error. I'll try the raster brick approach and see if it helps but i still think it could be something to do with pixel overlap - landsat scenes are not square but the extents are 'rectangular' so there are many no data pixels in a scene that could cause issues. I'll try a few more things ... – Leah Wasser Jun 01 '17 at 17:45
  • nope! as i suspected it is not related to being a brick vs stack. i still think it is something to do with pixel overlap / missing data . – Leah Wasser Jun 01 '17 at 17:50
1

This is a bug. Cover with two multi-layered Raster* objects cannot write to disk. This can be seen in the simple example by setting 

rasterOptions(todisk=TRUE)

I have fixed this in version 2-6.1 (forthcoming)

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63