r-spatial / rgee

Google Earth Engine for R
https://r-spatial.github.io/rgee/
Other
677 stars 146 forks source link

Reducing image collection to image returns only one value per band #214

Closed DavidDHofmann closed 2 years ago

DavidDHofmann commented 2 years ago

At submit an issue, please attached the following information of your rgee session:

library(rgee)

# Initialize the Earth Engine module.
ee_Initialize()

# Print metadata for a DEM dataset.
print(ee$Image('USGS/SRTMGL1_003')$getInfo())

Description

Hi there!

I'm trying to download Sentinel-2 data for a specific date. For the date and area of interest, Earth Engine returns an image-collection with two images that overlap in extent. I then try to reduce the collection into a single image using a median-reducer. Ultimately, I clip the reduced image to my area of interest and download the data. Oddly enough, ee_print(reduced) shows that the image only contains one pixel per band. Nevertheless, the image looks fine when displayed on the map. Only when downloaded, the bands contain only one value per band. To me, it seems that the reducer does not apply the median per pixel across bands, but across the entire band. The coordinates of the final image also do not exactly match my specified aoi and I get a warning after download saying:

Warning
In .rasterFromGDAL(x, band = band, objecttype, ...) :
  data seems flipped. Consider using: flip(x, direction='y')

I'm not quite sure what I'm doing wrong. Any suggestions?

What I Did

# Load required packages
library(rgee)
library(terra)

# Init
ee_Initialize()

# Specify an area of interest
aoi <- ee$Geometry$Polygon(
  list(
      c(23.4976149678726, -19.649075296214)
    , c(23.4976149678726, -19.2680589256603)
    , c(23.93394016641, -19.2680589256603)
    , c(23.93394016641, -19.649075296214)
  )
)

# Query sentinel 2 data
collection <- ee$ImageCollection("COPERNICUS/S2_SR")$
  filterDate("2021-04-05", "2021-04-06")$
  filterBounds(aoi)$
  select("B2")

# Check it
ee_print(collection)
────────────────────────────────────────────────────────────── Earth Engine ImageCollection ──
ImageCollection Metadata:
 - Class                      : ee$ImageCollection
 - Number of Images           : 2
 - Number of Properties       : 23
 - Number of Pixels*          : 241142760
 - Approximate size*          : 3.36 GB
Image Metadata (img_index = 0):
 - ID                         : COPERNICUS/S2_SR/20210405T081559_20210405T083938_T34KGD
 - system:time_start          : 2021-04-05 08:45:41
 - system:time_end            : 2021-04-05 08:45:41
 - Number of Bands            : 1
 - Bands names                : B2
 - Number of Properties       : 81
 - Number of Pixels*          : 120571380
 - Approximate size*          : 1.68 GB
Band Metadata (img_band = 'B2'):
 - EPSG (SRID)                : WGS 84 / UTM zone 34S (EPSG:32734)
 - proj4string                : +proj=utm +zone=34 +south +datum=WGS84 +units=m +no_defs
 - Geotransform               : 10 0 699960 0 -10 7900000
 - Nominal scale (meters)     : 10
 - Dimensions                 : 10981 10980
 - Number of Pixels           : 120571380
 - Data type                  : INT
 - Approximate size           : 1.68 GB
 ──────────────────────────────────────────────────────────────────────────────────────────────
 NOTE: (*) Properties calculated considering a constant  geotransform and data type.
# There are two images on this day, let's take the median of the matching bands
reduced <- collection$
  mean()$
  clip(aoi)

# Check again
ee_print(reduced)
# Displayed maps look correct
Map$centerObject(reduced, zoom = 10)
Map$addLayer(reduced$select("B2"), list(max = 1000, min = 0), name = "B2")

image

# Downloaded data shows only one value in the band
ee_as_raster(reduced, dsn = "reduced.tif")
r <- rast("reduced.tif")
show(r)
dimensions  : 1, 1, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 23, 24, -20, -19  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : reduced.tif 
name        : B2 
plot(r)

image

csaybar commented 2 years ago

Hi @DavidDHofmann you have to define the scale argument:

reduced2 <- ee_as_raster(reduced, dsn = "reduced.tif", scale = 1000)

The warnings appear because GEE retrieved your query using positive y spacing and raster does not automatically support it. I recommend you use ee_as_stars to avoid this problem. https://github.com/evetion/GeoArrays.jl/issues/53

Warning
In .rasterFromGDAL(x, band = band, objecttype, ...) :
  data seems flipped. Consider using: flip(x, direction='y')
DavidDHofmann commented 2 years ago

Hi @csaybar Thank you for the swift reply! The code works when providing the scale argument. Does this mean that it is not possible to download the data at its native resolution?

csaybar commented 2 years ago

Hi @DavidDHofmann :)

In fact, you can!. The issue is that after applying a reducer collection$mean()$clip(aoi) you will lose the image geotransform properties.

ee_print(collection)

image

ee_print(collection$mean())

image

It is a consequence of GEE lazy evaluation. The client library API does not know if after applying your "selected reducer" (mean) you are conserving the image geometric properties it's something that will be evaluated at the server-side. My advice would be to store the properties before the reducer.

# Load required packages
library(rgee)
library(terra)

# Init
ee_Initialize()

# Specify an area of interest
aoi <- ee$Geometry$Polygon(
  list(
    c(23.4976149678726, -19.649075296214)
    , c(23.4976149678726, -19.2680589256603)
    , c(23.93394016641, -19.2680589256603)
    , c(23.93394016641, -19.649075296214)
  )
)

# Query sentinel 2 data
collection <- ee$ImageCollection("COPERNICUS/S2_SR")$
  filterDate("2021-04-05", "2021-04-06")$
  filterBounds(aoi)$
  select("B2")

collection$first()$projection()$nominalScale()$getInfo()
#10
DavidDHofmann commented 2 years ago

Ah okey cool, I wasn't aware of this. This also explains why the estimated image size becomes unreasonably small after applying the reudcer. Thank you very much for the clarification. The download works perfectly fine now :)