1

I am trying to extract a number of variables at a series of points (lat, lon, time) by intersecting multiple images and image collections with those points (e.g. Landsat images and elevation raster). The code below gives me Landsat pixel values at the points, but I don't know how to combine these values with those extracted from another raster (e.g. elev = ee.Image('USGS/NED')). My end goal is to get a table whose rows are my points and columns are Landsat band values, elevation value at that point, etc., but do that in a single query because I have quite a few rasters to extract from.

Here is my code:

# Sample list of query points
pointList = [
  ee.Feature(ee.Geometry.Point([-121.771399, 40.330908],'EPSG:4269'), {'point_id': 'p1', 'point_date': '2019-08-27'}),
  ee.Feature(ee.Geometry.Point([ -120.482291,39.245937],'EPSG:4269'), {'point_id': 'p2', 'point_date': '2013-8-19'}),
  ee.Feature(ee.Geometry.Point([-120.085182,37.566734],'EPSG:4269'), {'point_id': 'p3', 'point_date': '2014-04-25'})
]

points_fc = ee.FeatureCollection(pointList)


def maskL8sr(image):
  # Bits 3 and 5 are cloud shadow and cloud, respectively.
  cloudShadowBitMask = (1 << 3)
  cloudsBitMask = (1 << 5)
  # Get the pixel QA band.
  qa = image.select('pixel_qa')
  # Both flags should be set to zero, indicating clear conditions.
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \
                 .And(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return image.updateMask(mask)


def extract(point):    
  ref_date = ee.Date(point.get('point_date'))
  start_date =  ref_date.advance(-1, 'month')
  end_date =  ref_date.advance(1, 'month')

  # Extract image collection at the desired point and time period
  L8_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
                    .filterBounds(point.geometry()) \
                    .filterDate(start_date, end_date) \
                    .map(maskL8sr)
  projection = L8_col.first().projection() # projection of first image
  median = L8_col.median()
    
  # Wrap the single image in an ImageCollection for mapping
  collection = ee.ImageCollection([median]); 
  

  def reduction(img):
    pix_val = img.reduceRegion(
              reducer= ee.Reducer.first(),
              geometry= point.geometry(),
              crs= projection,
              scale= projection.nominalScale(),
            )
    return point.set(pix_val)
  
  # Map over the collection and extract pixel values
  extracted = collection.map(reduction)
  return ee.FeatureCollection(extracted)

# Flatten the result
points_FC = points_fc.map(extract).flatten()

print(points_FC.getInfo())

1 Answer 1

2

When you want data from multiple images, a straightforward approach is to use addBands to combine the two images into one.

ned = ee.Image('USGS/NED')
pix_val = img.addBands(ned).reduceRegion(...

In some cases, you might need to rename the bands of one of the images, but that is not necessary here.


This isn't related to your question, but I notice you're going about this in a roundabout way:

  # Wrap the single image in an ImageCollection for mapping
  collection = ee.ImageCollection([median]); 
  

  def reduction(img):
    pix_val = img.reduceRegion(
              reducer= ee.Reducer.first(),
              geometry= point.geometry(),
              crs= projection,
              scale= projection.nominalScale(),
            )
    return point.set(pix_val)
  
  # Map over the collection and extract pixel values
  extracted = collection.map(reduction)
  return ee.FeatureCollection(extracted)

Unless you're intending to eventually generate multiple rows, there is no need to create a collection and call map here at all; just directly return the feature:

  pix_val = img.addBands(ned).reduceRegion(
            reducer= ee.Reducer.first(),
            geometry= point.geometry(),
            crs= projection,
            scale= projection.nominalScale(),
          )
  return point.set(pix_val)

With this change to extract to return a feature instead of a collection, you'll correspondingly need to remove the .flatten() after .map(extract).

4
  • This works great. Will addBands work if I want to extract all pixels in a neighborhood of the point? From your previous answer, will it be problematic if the added images have a different scale and projection? Basically, I add the additional bands, then use img_added.neighborhoodToArray(ee.Kernel.square(1), 999); and then reduce the neighborhood to extract all pixel values within 3x3 pixels (at the resolution of Landsat median, i.e. 90mx90m) around the point. Will this work if I use the scale and projection from the first Landsat image used to calculate the median within reduceRegion?
    – Mohamad
    Jul 22, 2021 at 14:57
  • 1
    "Will addBands work if I want to extract all pixels in a neighborhood of the point?": The neighbor will be defined in terms of one particular (re)projection as specified in reduceRegion; you'd have to do them separately and combine the resulting dictionaries if you want the neighborhoods of original pixels in two images (but notice that these may be rotated as well as at different scales).
    – Kevin Reid
    Jul 22, 2021 at 15:04
  • So addBands works great if I want just the value at the point from multiple images, but if I want more than one pixel using neighborhoodToArray, I shouldn't use addBands and instead, do separate reduceRegion extractions, and even then, the pixel grids extracted from the different images will not exactly align with each other (may be rotated and different scales). Is that a correct summary?
    – Mohamad
    Jul 22, 2021 at 15:22
  • 1
    @Mohamad Yes, except that "shouldn't" depends on what your goal is in getting the neighborhood. If you want data on the same grid, then a single reduceRegion is appropriate, and the problem is in choosing an appropriate grid given the limited resolution of the input images. If you want original image pixels, you have to deal with the fact that the images you're starting with don't have the same grid.
    – Kevin Reid
    Jul 22, 2021 at 18:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.