Sampling areas based on water mask using Google Earth Engine

Question

I am trying to sample a large image and I only want to keep points that fall within a water mask.

I can successfully sample small numbers of points by using a mask, but the algorithm chokes at larger numbers of pixels. This is surprising because I specify the numPixels argument as a large number, but the resultant featurecollection only consists of a small number of points.

Here is my .sample() call with the masking operation:

var water_mask = ee.Image("JRC/GSW1_3/GlobalSurfaceWater")
                    .select(['max_extent'],['is_water'])

var montana = ee.FeatureCollection("TIGER/2018/States")
                   .filter(ee.Filter.eq("NAME","Montana"))

var my_img = img1.addBands(img2).mask(water_mask)


// Call to .sample()
var samples = my_img.sample({'region':montana,
                          'projection':"EPSG:4326",
                          'scale':30,
                          'numPixels':500})
print(samples.size())
// Even though numPixels is specified as 500, the size of 'samples' is only 7

This behavior is synonymous to sampling an unmasked image and filtering the resultant featurecollection by the mask.

For example,

// First add the water_mask bands to the image, so that they can be sampled
var synonymous = my_img.addBands(water_mask)
                 .addBands(ndvi)

                 // Sample masked image - note numPixels
                 .sample({'region':montana,
                          'projection':"EPSG:4326",
                          'scale':30,
                          'numPixels':500})
                          
                 // Filter feature collection by condition
                 .filter(ee.Filter.eq('is_water',1))

print(synonymous.size())
// This results in only 7 points as well

Is there a more efficient way to do this operation?

Here is a link to my code: https://code.earthengine.google.com/195b11ed5981b695bfe9bbe2da41acac

Answer 1

The docs for sample() explicitly state:

Note that the default behavior is to drop features that intersect masked pixels, which result in null-valued properties.

Even though you've only got 1 class, what you're trying to do is a stratified sampling, so use stratifiedSample.

var samples = lst.addBands(water_mask)
                .addBands(ndvi)
                .updateMask(water_mask)
                .stratifiedSample({
                  numPoints: 500,
                  classBand: 'is_water',
                  region: montana,
                  projection: "EPSG:4326",
                  scale :30
                })

Note that Montana at 30m is pretty big (422 million pixels), so you might have to do this with an export.

Also, don't clip; you're already using the geometry for a region, so clipping the images wastes time/memory.