stars analogous of raster::aggregate

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

stars analogous of raster::aggregate

Julian M. Burgos
Dear all,

The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10:

library(raster)

rst1 <-  raster(system.file("tif/L7_ETMs.tif", package = "stars"))

rst2 <- aggregate(rst1, fact = 10, fun = mean)

> res(rst1)
[1] 28.5 28.5

> res(rst2)
[1] 285 285

I am trying to do the same thing with a stars object.  The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons.  I could do something like this:

library(stars)
st1 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>%
  slice(band, 1)
grid <- st_make_grid(st_bbox(st1) %>%
                     st_as_sfc(),
                     cellsize = 285)

st2 <- aggregate(st1, grid, mean)

But the resulting st2 object is not a raster (i.e. does not X and Y dimensions), but instead has a single dimension and contains the geometries of the grid polygons.

> st1
stars object with 2 dimensions and 1 attribute
attribute(s):
  L7_ETMs.tif
 Min.   : 47.00
 1st Qu.: 67.00
 Median : 78.00
 Mean   : 79.15
 3rd Qu.: 89.00
 Max.   :255.00
dimension(s):
  from  to  offset delta                       refsys point values
x    1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
y    1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
> st2
stars object with 1 dimensions and 1 attribute
attribute(s):
  L7_ETMs.tif
 Min.   : 57.76
 1st Qu.: 69.22
 Median : 79.04
 Mean   : 79.09
 3rd Qu.: 87.63
 Max.   :123.51
dimension(s):
         from   to offset delta                       refsys point
geometry    1 1260     NA    NA UTM Zone 25, Southern Hem... FALSE
                                                                    values
geometry POLYGON ((288776.3 9110729,...,...,POLYGON ((298466.3 9120704,...
>

So I am stuck, wondering if perhaps st_apply is the answer.  What would be the best way to replicate what the raster::aggregate function does using stars?  Any ideas?

Many thanks,

Julian

--
Julian Mariano Burgos, PhD
Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/
Marine and Freshwater Research Institute
Botnsjávarsviðs / Demersal Division
  Fornubúðir 5, IS-220 Hafnarfjörður, Iceland
www.hafogvatn.is
Sími/Telephone : +354-5752037
Netfang/Email: [hidden email]
rr

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: stars analogous of raster::aggregate

Micha Silver


On 08/09/2020 19:33, Julian M. Burgos wrote:
Dear all,

The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10:

library(raster)

rst1 <-  raster(system.file("tif/L7_ETMs.tif", package = "stars"))

rst2 <- aggregate(rst1, fact = 10, fun = mean)

res(rst1)
[1] 28.5 28.5

res(rst2)
[1] 285 285

I am trying to do the same thing with a stars object.  The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons.  I could do something like this:

Probably st_warp() is what you want:


l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
l7 = read_stars(l7_file)

l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs = st_crs(l7))


stars::st_dimensions(l7)
     from  to  offset delta                       refsys point values   
x       1 349  288776  28.5 +proj=utm +zone=25 +south... FALSE   NULL [x]
y       1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE   NULL [y]
band    1   6      NA    NA                           NA    NA   NULL   


stars::st_dimensions(l7_lowres)
     from  to  offset delta                       refsys point values   
x       1 111  288776    90 +proj=utm +zone=25 +south...    NA   NULL [x]
y       1 112 9120761   -90 +proj=utm +zone=25 +south...    NA   NULL [y]
band    1   6      NA    NA                           NA    NA   NULL  


-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://orcid.org/0000-0002-1128-1325

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: stars analogous of raster::aggregate

edzer


On 9/8/20 8:47 PM, Micha Silver wrote:

>
> On 08/09/2020 19:33, Julian M. Burgos wrote:
>> Dear all,
>>
>> The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10:
>>
>> library(raster)
>>
>> rst1 <-  raster(system.file("tif/L7_ETMs.tif", package = "stars"))
>>
>> rst2 <- aggregate(rst1, fact = 10, fun = mean)
>>
>>> res(rst1)
>> [1] 28.5 28.5
>>
>>> res(rst2)
>> [1] 285 285
>>
>> I am trying to do the same thing with a stars object.  The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons.  I could do something like this:
>
> Probably st_warp() is what you want:
>
>
> l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
> l7 = read_stars(l7_file)
>
> l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs = st_crs(l7))
>
>
> stars::st_dimensions(l7)
>       from  to  offset delta                       refsys point values
> x       1 349  288776  28.5 +proj=utm +zone=25 +south... FALSE NULL [x]
> y       1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL [y]
> band    1   6      NA    NA                           NA    NA NULL
>
>
> stars::st_dimensions(l7_lowres)
>       from  to  offset delta                       refsys point values
> x       1 111  288776    90 +proj=utm +zone=25 +south...    NA NULL [x]
> y       1 112 9120761   -90 +proj=utm +zone=25 +south...    NA NULL [y]
> band    1   6      NA    NA                           NA    NA NULL

Yes; if you want to get similar behaviour as raster, choose a cell size
that is an exact multiple of the origin's cell size, use_gdal = TRUE,
and method = "average"; this currently seems to only work for single
band rasters; along these lines:

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
l7 = read_stars(l7_file)
(dest = st_as_stars(st_bbox(l7), dx = 90, dy = 90))
# stars object with 2 dimensions and 1 attribute
# attribute(s):
#     values
#  Min.   :0
#  1st Qu.:0
#  Median :0
#  Mean   :0
#  3rd Qu.:0
#  Max.   :0
# dimension(s):
#   from  to  offset delta                       refsys point values x/y
# x    1 111  288776    90 UTM Zone 25, Southern Hem...    NA   NULL [x]
# y    1 112 9120761   -90 UTM Zone 25, Southern Hem...    NA   NULL [y]
(l7_lowres = st_warp(src = l7[,,,1], dest = dest, use_gdal = TRUE,
method = "average"))
# stars object with 2 dimensions and 1 attribute
# attribute(s):
#  file14bb5c7321dc.tif
#  Min.   : 56.81
#  1st Qu.: 68.75
#  Median : 79.00
#  Mean   : 79.25
#  3rd Qu.: 88.62
#  Max.   :207.44
# dimension(s):
#   from  to  offset delta                       refsys point values x/y
# x    1 111  288776    90 UTM Zone 25, Southern Hem... FALSE   NULL [x]
# y    1 112 9120761   -90 UTM Zone 25, Southern Hem... FALSE   NULL [y]

>
>
> --
> Micha Silver
> Ben Gurion Univ.
> Sde Boker, Remote Sensing Lab
> cell: +972-523-665918
> https://orcid.org/0000-0002-1128-1325
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48149 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: stars analogous of raster::aggregate

Julian M. Burgos
Thanks Edzer and Micha, this was very helpful.  st_warp with use_gdal=TRUE worked like a charm.

My best,

Julian


Edzer Pebesma writes:

> On 9/8/20 8:47 PM, Micha Silver wrote:
>> On 08/09/2020 19:33, Julian M. Burgos wrote:
>>> Dear all,
>>>
>>> The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10:
>>>
>>> library(raster)
>>>
>>> rst1 <-  raster(system.file("tif/L7_ETMs.tif", package = "stars"))
>>>
>>> rst2 <- aggregate(rst1, fact = 10, fun = mean)
>>>
>>>> res(rst1)
>>> [1] 28.5 28.5
>>>
>>>> res(rst2)
>>> [1] 285 285
>>>
>>> I am trying to do the same thing with a stars object.  The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons.  I could do something like this:
>> Probably st_warp() is what you want:
>>
>> l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
>> l7 = read_stars(l7_file)
>> l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs =
>> st_crs(l7))
>>
>> stars::st_dimensions(l7)
>>       from  to  offset delta                       refsys point values
>> x       1 349  288776  28.5 +proj=utm +zone=25 +south... FALSE NULL [x]
>> y       1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL [y]
>> band    1   6      NA    NA                           NA    NA NULL
>>
>> stars::st_dimensions(l7_lowres)
>>       from  to  offset delta                       refsys point values
>> x       1 111  288776    90 +proj=utm +zone=25 +south...    NA NULL [x]
>> y       1 112 9120761   -90 +proj=utm +zone=25 +south...    NA NULL [y]
>> band    1   6      NA    NA                           NA    NA NULL
>
> Yes; if you want to get similar behaviour as raster, choose a cell
> size that is an exact multiple of the origin's cell size, use_gdal =
> TRUE, and method = "average"; this currently seems to only work for
> single band rasters; along these lines:
>
> library(stars)
> # Loading required package: abind
> # Loading required package: sf
> # Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
> l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
> l7 = read_stars(l7_file)
> (dest = st_as_stars(st_bbox(l7), dx = 90, dy = 90))
> # stars object with 2 dimensions and 1 attribute
> # attribute(s):
> #     values
> #  Min.   :0
> #  1st Qu.:0
> #  Median :0
> #  Mean   :0
> #  3rd Qu.:0
> #  Max.   :0
> # dimension(s):
> #   from  to  offset delta                       refsys point values x/y
> # x    1 111  288776    90 UTM Zone 25, Southern Hem...    NA   NULL [x]
> # y    1 112 9120761   -90 UTM Zone 25, Southern Hem...    NA   NULL [y]
> (l7_lowres = st_warp(src = l7[,,,1], dest = dest, use_gdal = TRUE,
> method = "average"))
> # stars object with 2 dimensions and 1 attribute
> # attribute(s):
> #  file14bb5c7321dc.tif
> #  Min.   : 56.81
> #  1st Qu.: 68.75
> #  Median : 79.00
> #  Mean   : 79.25
> #  3rd Qu.: 88.62
> #  Max.   :207.44
> # dimension(s):
> #   from  to  offset delta                       refsys point values x/y
> # x    1 111  288776    90 UTM Zone 25, Southern Hem... FALSE   NULL [x]
> # y    1 112 9120761   -90 UTM Zone 25, Southern Hem... FALSE   NULL [y]
>
>>
>> -- Micha Silver
>> Ben Gurion Univ.
>> Sde Boker, Remote Sensing Lab
>> cell: +972-523-665918
>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F0000-0002-1128-1325&amp;data=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294&amp;sdata=K5kMOHf6dka1eNJrVRZsgDe0UmxTEhBRU40C%2FPpRXVE%3D&amp;reserved=0
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294&amp;sdata=KwymhTSAHdcBDVffZHhp2xZ%2BQUzDKh%2B%2BN2DP1yLtGY4%3D&amp;reserved=0
>>


--
Julian Mariano Burgos, PhD
Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/
Marine and Freshwater Research Institute
Botnsjávarsviðs / Demersal Division
  Fornubúðir 5, IS-220 Hafnarfjörður, Iceland
www.hafogvatn.is
Sími/Telephone : +354-5752037
Netfang/Email: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: stars analogous of raster::aggregate

edzer
It now also works for multi-band rasters (fixed in sf when installed
from github): https://github.com/r-spatial/stars/issues/320

On 9/9/20 1:10 PM, Julian M. Burgos wrote:

> Thanks Edzer and Micha, this was very helpful.  st_warp with use_gdal=TRUE worked like a charm.
>
> My best,
>
> Julian
>
>
> Edzer Pebesma writes:
>
>> On 9/8/20 8:47 PM, Micha Silver wrote:
>>> On 08/09/2020 19:33, Julian M. Burgos wrote:
>>>> Dear all,
>>>>
>>>> The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10:
>>>>
>>>> library(raster)
>>>>
>>>> rst1 <-  raster(system.file("tif/L7_ETMs.tif", package = "stars"))
>>>>
>>>> rst2 <- aggregate(rst1, fact = 10, fun = mean)
>>>>
>>>>> res(rst1)
>>>> [1] 28.5 28.5
>>>>
>>>>> res(rst2)
>>>> [1] 285 285
>>>>
>>>> I am trying to do the same thing with a stars object.  The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons.  I could do something like this:
>>> Probably st_warp() is what you want:
>>>
>>> l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
>>> l7 = read_stars(l7_file)
>>> l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs =
>>> st_crs(l7))
>>>
>>> stars::st_dimensions(l7)
>>>        from  to  offset delta                       refsys point values
>>> x       1 349  288776  28.5 +proj=utm +zone=25 +south... FALSE NULL [x]
>>> y       1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL [y]
>>> band    1   6      NA    NA                           NA    NA NULL
>>>
>>> stars::st_dimensions(l7_lowres)
>>>        from  to  offset delta                       refsys point values
>>> x       1 111  288776    90 +proj=utm +zone=25 +south...    NA NULL [x]
>>> y       1 112 9120761   -90 +proj=utm +zone=25 +south...    NA NULL [y]
>>> band    1   6      NA    NA                           NA    NA NULL
>>
>> Yes; if you want to get similar behaviour as raster, choose a cell
>> size that is an exact multiple of the origin's cell size, use_gdal =
>> TRUE, and method = "average"; this currently seems to only work for
>> single band rasters; along these lines:
>>
>> library(stars)
>> # Loading required package: abind
>> # Loading required package: sf
>> # Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
>> l7_file = system.file("tif/L7_ETMs.tif", package = "stars")
>> l7 = read_stars(l7_file)
>> (dest = st_as_stars(st_bbox(l7), dx = 90, dy = 90))
>> # stars object with 2 dimensions and 1 attribute
>> # attribute(s):
>> #     values
>> #  Min.   :0
>> #  1st Qu.:0
>> #  Median :0
>> #  Mean   :0
>> #  3rd Qu.:0
>> #  Max.   :0
>> # dimension(s):
>> #   from  to  offset delta                       refsys point values x/y
>> # x    1 111  288776    90 UTM Zone 25, Southern Hem...    NA   NULL [x]
>> # y    1 112 9120761   -90 UTM Zone 25, Southern Hem...    NA   NULL [y]
>> (l7_lowres = st_warp(src = l7[,,,1], dest = dest, use_gdal = TRUE,
>> method = "average"))
>> # stars object with 2 dimensions and 1 attribute
>> # attribute(s):
>> #  file14bb5c7321dc.tif
>> #  Min.   : 56.81
>> #  1st Qu.: 68.75
>> #  Median : 79.00
>> #  Mean   : 79.25
>> #  3rd Qu.: 88.62
>> #  Max.   :207.44
>> # dimension(s):
>> #   from  to  offset delta                       refsys point values x/y
>> # x    1 111  288776    90 UTM Zone 25, Southern Hem... FALSE   NULL [x]
>> # y    1 112 9120761   -90 UTM Zone 25, Southern Hem... FALSE   NULL [y]
>>
>>>
>>> -- Micha Silver
>>> Ben Gurion Univ.
>>> Sde Boker, Remote Sensing Lab
>>> cell: +972-523-665918
>>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F0000-0002-1128-1325&amp;data=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294&amp;sdata=K5kMOHf6dka1eNJrVRZsgDe0UmxTEhBRU40C%2FPpRXVE%3D&amp;reserved=0
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294&amp;sdata=KwymhTSAHdcBDVffZHhp2xZ%2BQUzDKh%2B%2BN2DP1yLtGY4%3D&amp;reserved=0
>>>
>
>
> --
> Julian Mariano Burgos, PhD
> Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/
> Marine and Freshwater Research Institute
> Botnsjávarsviðs / Demersal Division
>    Fornubúðir 5, IS-220 Hafnarfjörður, Iceland
> www.hafogvatn.is
> Sími/Telephone : +354-5752037
> Netfang/Email: [hidden email]
>

--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48149 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo