Long extract time for multiband raster

SFrav commented 3 years ago

There's probably something basic that I'm doing wrong here, but I have two raster stacks that are derived from the same primary source (same resolution, extent and number of layers). With a vector layer of ~40k features, one stack takes minutes and the other stack took about 18 hours. The only difference that I can see is that the one that took way longer has two bands.

I imported that one like so: st_laiSNAP <- stack(list.files(path = 'LAI_S2_SNAP/', pattern = "\\.tif$", full.names = T), bands = 1)

The output of str on the problem stack is below

The differences I see are:

@ blockrows/cols : int 10980 (the quicker stack is int 256)
@ nbands : int 2 (the quicker stack has 1)

The exactextractr output of the problem stack is just as expected. No problems there.

With the same resolution, extent, n layers and features, why is there such a big difference in processing time?

Exactextractr version: 0.2.1

Problem stack: Formal class 'RasterStack' [package "raster"] with 11 slots ..@ filename: chr "" ..@ layers :List of 6 .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A006069_20180505T112304_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi TRUE .. .. .. .. .. ..@ min : num -0.672 .. .. .. .. .. ..@ max : num 13.9 .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006069_20180505T112304_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A006355_20180525T112111_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006355_20180525T112111_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A006398_20180528T113438_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006398_20180528T113438_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A006541_20180607T113457_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006541_20180607T113457_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A006927_20180704T112112_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006927_20180704T112112_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "\L1C_T30UWG_A015764_20180629T112537_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A015764_20180629T112537_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() ..@ title : chr(0) ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. ..@ xmin: num 5e+05 .. .. ..@ xmax: num 609780 .. .. ..@ ymin: num 6090240 .. .. ..@ ymax: num 6200040 ..@ rotated : logi FALSE ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. ..@ geotrans: num(0) .. .. ..@ transfun:function ()
..@ ncols : int 10980 ..@ nrows : int 10980 ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" ..@ history : list() ..@ z : list()

Problem stack Formal class 'RasterStack' [package "raster"] with 11 slots ..@ filename: chr "" ..@ layers :List of 6 .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A006069_20180505T112304_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi TRUE .. .. .. .. .. ..@ min : num -0.672 .. .. .. .. .. ..@ max : num 13.9 .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006069_20180505T112304_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A006355_20180525T112111_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006355_20180525T112111_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A006398_20180528T113438_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006398_20180528T113438_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A006541_20180607T113457_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006541_20180607T113457_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A006927_20180704T112112_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A006927_20180704T112112_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() .. ..$ :Formal class 'RasterLayer' [package "raster"] with 12 slots .. .. .. ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots .. .. .. .. .. ..@ name : chr "C:\Users\Simon\Documents\GAAFS\Ideas\EOcropModel\LAI_S2_SNAP\L1C_T30UWG_A015764_20180629T112537_10m_lai.tif" .. .. .. .. .. ..@ datanotation: chr "FLT4S" .. .. .. .. .. ..@ byteorder : chr "little" .. .. .. .. .. ..@ nodatavalue : num -Inf .. .. .. .. .. ..@ NAchanged : logi FALSE .. .. .. .. .. ..@ nbands : int 2 .. .. .. .. .. ..@ bandorder : chr "BIL" .. .. .. .. .. ..@ offset : int 0 .. .. .. .. .. ..@ toptobottom : logi TRUE .. .. .. .. .. ..@ blockrows : int 10980 .. .. .. .. .. ..@ blockcols : int 10980 .. .. .. .. .. ..@ driver : chr "gdal" .. .. .. .. .. ..@ open : logi FALSE .. .. .. ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ offset : num 0 .. .. .. .. .. ..@ gain : num 1 .. .. .. .. .. ..@ inmemory : logi FALSE .. .. .. .. .. ..@ fromdisk : logi TRUE .. .. .. .. .. ..@ isfactor : logi FALSE .. .. .. .. .. ..@ attributes: list() .. .. .. .. .. ..@ haveminmax: logi FALSE .. .. .. .. .. ..@ min : num Inf .. .. .. .. .. ..@ max : num -Inf .. .. .. .. .. ..@ band : int 1 .. .. .. .. .. ..@ unit : chr "" .. .. .. .. .. ..@ names : chr "L1C_T30UWG_A015764_20180629T112537_10m_lai" .. .. .. ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots .. .. .. .. .. ..@ type : chr(0) .. .. .. .. .. ..@ values : logi(0) .. .. .. .. .. ..@ color : logi(0) .. .. .. .. .. ..@ names : logi(0) .. .. .. .. .. ..@ colortable: logi(0) .. .. .. ..@ title : chr(0) .. .. .. ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. .. .. .. ..@ xmin: num 5e+05 .. .. .. .. .. ..@ xmax: num 609780 .. .. .. .. .. ..@ ymin: num 6090240 .. .. .. .. .. ..@ ymax: num 6200040 .. .. .. ..@ rotated : logi FALSE .. .. .. ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. .. .. .. ..@ geotrans: num(0) .. .. .. .. .. ..@ transfun:function ()
.. .. .. ..@ ncols : int 10980 .. .. .. ..@ nrows : int 10980 .. .. .. ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. .. .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" .. .. .. ..@ history : list() .. .. .. ..@ z : list() ..@ title : chr(0) ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots .. .. ..@ xmin: num 5e+05 .. .. ..@ xmax: num 609780 .. .. ..@ ymin: num 6090240 .. .. ..@ ymax: num 6200040 ..@ rotated : logi FALSE ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots .. .. ..@ geotrans: num(0) .. .. ..@ transfun:function ()
..@ ncols : int 10980 ..@ nrows : int 10980 ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot .. .. ..@ projargs: chr "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" ..@ history : list() ..@ z : list()

dbaston commented 3 years ago

I'd first suggest updating to the most current version. From the changelog, it looks like this may have been addressed in June with version 0.4.0. If the problem persists, can you please share enough data and code that I can reproduce it?

SFrav commented 3 years ago

Thanks for the quick reply. I updated, but now I get the following error on running for any stack

Error in sf::st_geometry_type(y, by_geometry = FALSE) : unused argument (by_geometry = FALSE)

I'll reboot and try again. From what I can see, there are no changes in syntax that would have made a difference.

SFrav commented 3 years ago

This error occurs with the basic example as well (below). I probably have other packages to update to get this to work.

I'll check dependency versions and try again. Most likely in a few weeks.

rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
exact_extract(rast, poly, 'mean')

dbaston commented 3 years ago

Looks like you need sf version 0.9.0 or greater. I've just specified this in DESCRIPTION so install.packages will take care of it by itself in the future.

dbaston commented 3 years ago

Closing for now; please re-open if you're seeing the performance problem with a current version.

SFrav commented 3 years ago

This is still an issue with the current version. I could send an 8mb subsection of the raster and the vector layer that I'm working with. Would that help?

dbaston commented 3 years ago

Yes please. Whatever data/code is needed to reproduce would be helpful.

dbaston commented 3 years ago

Are you able to provide anything to reproduce this, @SFrav ?

SFrav commented 3 years ago

@dbaston, all the reproducible examples that I've produced resolve the processing time issue. Cropping the problem layers seems to resolve it - even with a large buffer around them. The original files are 900mb+

Even reducing the extent by a small amount resolves it. I'll close for now.

dbaston commented 3 years ago

Thanks for trying. I'm fine with debugging this on the original layers, but I understand you may not be able to share them.

SFrav commented 3 years ago

It's all open data, so no worries there. Coming your way in a few hrs.

dbaston commented 3 years ago

It looks like the issue is somewhat particular to this dataset. If you look at the L1C_T30UWG_A006069_20180505T112304_10m_lai.tif with rgdal::GDALinfo, you can see that it has a block size of 10980x10980 pixels (the entire image!). That means that whenever we read from the raster, GDAL will read in 120 million pixels, even if we only need the value of one. When we're reading from a single layer, that's not a huge problem, because GDAL will cache the values of those pixels, so subsequent reads will be fast:

> lyr <- raster('L1C_T30UWG_A006069_20180505T112304_10m_lai.tif')
> lyr <- readStart(lyr)
> replicate(3, system.time({
+   getValuesBlock(lyr, row = 10, nrows=10, col = 10, ncols = 10)
+ }))
            [,1] [,2]  [,3]
user.self  0.158    0 0.000
sys.self   0.183    0 0.001
elapsed    0.341    0 0.001
user.child 0.000    0 0.000
sys.child  0.000    0 0.000

With a stack of N layers, every time we try to read a location we are pulling in 120 million * N pixels. This doesn't fit in the cache, so we end up reading from disk every time, and subsequent reads are not significantly faster than the initial one.

> stk <- stack(rep.int('L1C_T30UWG_A006069_20180505T112304_10m_lai.tif', 6))
> stk <- readStart(stk)
> replicate(3, system.time({
+   getValuesBlock(stk, row = 10, nrows=10, col = 10, ncols = 10)
+ }))
            [,1]  [,2]  [,3]
user.self  1.899 1.841 1.807
sys.self   1.349 0.841 0.851
elapsed    3.248 2.682 2.658
user.child 0.000 0.000 0.000
sys.child  0.000 0.000 0.000

The second dataset (T30UWG_20180505T112109_lai.tif) uses a much more typical block size of 256x256 pixels. So when we read from a stack of N layers, we are keeping only 65,536 * N pixels in cache. We can have a quite a few layers in the stack before the cache is overwhelmed - and even if it is overwhelmed, it doesn't take long to read 65k pixels from disk.

> lyr <- raster('T30UWG_20180505T112109_lai.tif')
> lyr <- readStart(lyr)
> replicate(3, system.time({
+   getValuesBlock(lyr, row = 10, nrows=10, col = 10, ncols = 10)
+ }))
            [,1]  [,2]  [,3]
user.self  0.000 0.000 0.001
sys.self   0.000 0.000 0.000
elapsed    0.001 0.001 0.001
user.child 0.000 0.000 0.000
sys.child  0.000 0.000 0.000
> 
> stk <- stack(rep.int('T30UWG_20180505T112109_lai.tif', 6))
> stk <- readStart(stk)
> replicate(3, system.time({
+   getValuesBlock(stk, row = 10, nrows=10, col = 10, ncols = 10)
+ }))
            [,1]  [,2]  [,3]
user.self  0.000 0.001 0.002
sys.self   0.006 0.001 0.000
elapsed    0.006 0.002 0.002
user.child 0.000 0.000 0.000
sys.child  0.000 0.000 0.000

exactextractr is making the assumption that if you feed it RasterStack, it's advantageous to read from every layer of the stack at once before moving on to the next polygon. This is a case where that assumption is bad, and it would be much faster to loop over each layer of the RasterStack individually.

SFrav commented 3 years ago

Ok. That was the other difference I mentioned in the first post. I don't quite understand how block size is set. From what you've said above I understand that it defines the number of pixels bought into memory by default - or something along those lines. Is it set when exporting using, say, gdal? Is there a way to set block size when importing in R or Python? If I can understand a bit more then I can post an issue on the ESA SNAP repository.

dbaston commented 3 years ago

The block size is set when writing the file. By default, GDAL will use a block size of a single row...in this case, 10980x1. This is what you'd get from gdal_translate L1C_T30UWG_A006069_20180505T112304_10m_lai.tif test.tif

You can customize that with, e.g.,

gdal_translate L1C_T30UWG_A006069_20180505T112304_10m_lai.tif test.tif -co TILED=YES

which gives 256x256, or

gdal_translate L1C_T30UWG_A006069_20180505T112304_10m_lai.tif test.tif -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512

which gives 512x512. It does seem odd to me if SNAP is writing these with a block size of 10980x10980.

One last thing. GDAL must read an entire block at a time if the GeoTIFF is compressed, but if it is uncompressed it can read only the requested window. I noticed that L1C_T30UWG_A006069_20180505T112304_10m_lai.tif is uncompressed, so if you set Sys.setenv(GTIFF_DIRECT_IO='ON') in R before opening the file, you can avoid reading the entire block and get much better performance.

isciences / exactextractr

Long extract time for multiband raster #45