Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Long extract time for multiband raster #45

Closed
SFrav opened this issue Dec 10, 2020 · 14 comments
Closed

Long extract time for multiband raster #45

SFrav opened this issue Dec 10, 2020 · 14 comments
Labels
awaiting feedback Further information is requested performance

Comments

@SFrav
Copy link

SFrav commented Dec 10, 2020

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
Copy link
Member

dbaston commented Dec 10, 2020

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
Copy link
Author

SFrav commented Dec 10, 2020

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
Copy link
Author

SFrav commented Dec 10, 2020

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 added a commit that referenced this issue Dec 11, 2020
Need 0.9.0 for by_geometry argument to st_geometry_type

References #45
@dbaston
Copy link
Member

dbaston commented Dec 11, 2020

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
Copy link
Member

dbaston commented Dec 15, 2020

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

@dbaston dbaston closed this as completed Dec 15, 2020
@SFrav
Copy link
Author

SFrav commented Dec 15, 2020

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
Copy link
Member

dbaston commented Dec 15, 2020

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

@dbaston dbaston reopened this Dec 15, 2020
@dbaston dbaston added awaiting feedback Further information is requested performance labels Jan 13, 2021
@dbaston
Copy link
Member

dbaston commented Mar 9, 2021

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

@SFrav
Copy link
Author

SFrav commented Mar 15, 2021

@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.

@SFrav SFrav closed this as completed Mar 15, 2021
@dbaston
Copy link
Member

dbaston commented Mar 15, 2021

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
Copy link
Author

SFrav commented Mar 15, 2021

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

@dbaston
Copy link
Member

dbaston commented Mar 16, 2021

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
Copy link
Author

SFrav commented Mar 16, 2021

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
Copy link
Member

dbaston commented Mar 16, 2021

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.

dbaston added a commit to dbaston/exactextract that referenced this issue May 16, 2024
Need 0.9.0 for by_geometry argument to st_geometry_type

References isciences/exactextractr#45
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
awaiting feedback Further information is requested performance
Projects
None yet
Development

No branches or pull requests

2 participants