Using rgee to create histograms from Google Earth Engine samples - by region and band

Question

I'm trying to use Google Earth Engine for the first time, via the package rgee. I want to randomly sample pixels within an image that has three bands, across a set of different regions. With those sampled data I wanto to create a plot containing histograms for each sample region, for a given band. So I will have three plots (one for each band).

In my case a region represents a specific type of habitat within the image, so I have a polygon for each habitat type. Here's my code:

# Load packages and initialize connection to GEE
library(rgee)
ee_Initialize()

###############################################################
# Create a geometry for each habitat type that I want to sample
###############################################################

unimproved_grazing = ee$Geometry$Polygon(
  list(
    c(-1.345864517292552, 59.97338078318311),
    c(-1.345864517292552, 59.97237144577356),
    c(-1.3438904114575911, 59.97237144577356),
    c(-1.3438904114575911, 59.97338078318311)
  )
)

improved_cut = ee$Geometry$Polygon(
  list(
    c(-1.329604032461742, 59.92455155105425),
    c(-1.329604032461742, 59.923922477208166),
    c(-1.3280268935609851, 59.923922477208166),
    c(-1.3280268935609851, 59.92455155105425)
  )
)

heath = ee$Geometry$Polygon(
  list(
    c(-1.351189448199126, 59.94125887709115),
    c(-1.351189448199126, 59.940151826756974),
    c(-1.3490865963314502, 59.940151826756974),
    c(-1.3490865963314502, 59.94125887709115)
  )
)

arable = ee$Geometry$Polygon(
  list(
    c(-1.3263107736431534, 59.969782903189916),
    c(-1.3263107736431534, 59.96911710113049),
    c(-1.325023313316005, 59.96911710113049),
    c(-1.325023313316005, 59.969782903189916)
  )
)

##################################################
# Create feature collection of all habitat regions
##################################################

habitats = ee$FeatureCollection(
  list(
    ee$Feature(unimproved_grazing, list(name = "unimproved")),
    ee$Feature(improved_cut,       list(name = "improved")),
    ee$Feature(heath,              list(name = "heath")),
    ee$Feature(arable,             list(name = "arable"))  
  )
)

#####################################
# Load the Sentinel-1 ImageCollection
#####################################

sentinel1 = ee$ImageCollection('COPERNICUS/S1_GRD')$
  # Filter date range - one year of images
  filterDate("2019-03-01", "2020-03-01")

# Filter by metadata properties
vvvh = sentinel1$
  filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VV'))$
  filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VH'))$
  filter(ee$Filter$eq('instrumentMode', 'IW'))$
  filter(ee$Filter$eq('resolution_meters', 10))

# Filter to get images from different look angles.
vhAscending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'ASCENDING'))
vhDescending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'DESCENDING'))

# Create a composite from means at different polarizations and look angles.
# Calculate mean of VH ascending polarisation
VH_Ascending_mean <- vhAscending$select("VH")$mean()
# Calculate mean of VH descending polarisation
VH_Descending_mean <- vhDescending$select("VH")$mean()
# Merge VV and VH and calculate mean
VV_Ascending_Descending_mean <- vhAscending$select("VV") %>%
  ee$ImageCollection$merge(vhDescending$select("VV")) %>%
  ee$ImageCollection$mean()

# Create single image containing bands of interest
collection <- ee$ImageCollection$fromImages(list(
  VH_Ascending_mean,
  VV_Ascending_Descending_mean,
  VH_Descending_mean
))$toBands()

######################################
# Sample pixels in each habitat region
######################################
mySamples = collection$sampleRegions(
                    collection = habitats,
                    scale= 10
                    )

I can export this to a csv, and get the samples for the three bands in my image. I'd like to have these samples classified by region. Also, I'd like to directly access the data without exporting so that I can create a set of histograms. My preference would to plot the histogram using ggplot2, along the lines of:

ggplot(mySamples, aes(sample)) +
  geom_histogram() +
  facet_grid(band)

Where each histogram covers multiple habitats, but for one image band. So there would be three plots in my example (VV, VH_1 and VH_2) each showing a histogram for the different sampling regions (habitats).

Answer 1

Taking advantage of the argument geometries=TRUE and ee_as_sf we can retrieve data with a single line of code. Also, it's a good practice to retrieve data in a batch way always it is possible.

######################################
# Sample pixels in each habitat region
######################################
ee_help(collection$sampleRegions)

for (index in 0:2) {
  habitat_f <- habitats %>% ee_get(index)
  mySamples <- collection$sampleRegions(
    collection = ee$Feature(habitat_f$first()),
    scale= 10,
    geometries=TRUE
  )
  cat(sprintf("Processing geom: %02d \n", index + 1))
  mySamples_sf <- ee_as_sf(mySamples)
  if (index == 0) {
    mySamples_sf_batch <- mySamples_sf   
  } else {
    mySamples_sf_batch <- rbind(mySamples_sf_batch, mySamples_sf)
  }
}

ggplot(mySamples_sf_batch, aes(X1_VV)) +
  geom_histogram() +
  facet_grid(~name)