1

I am trying to extract Landsat values at a number of points using Google Earth Engine (Python API). My points are given to me in the NAD83 system (as a CSV dataframe).

How can I extract a 3x3 grid of pixels (9 sets of spectral values) centered at the point from Landsat instead of just a single (30mx30m) pixel at the lat-lon point, while also enforcing the specific coordinate system (NAD83)? In other words, get all 30mx30m Landsat pixels within a 90mx90m area around my point?

Here is my code:

lat,lon = (40.330908, -121.771399) # coordinates in NAD83

point = ee.Geometry.Point([lon, lat])
landsat = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR') \
    .filterBounds(point) \
    .filterDate('2005-03-11', '2005-05-10')\
    .median()

extracted = landsat.reduceRegion(ee.Reducer.first(), point, scale = 30, crs ='EPSG:4269').getInfo()
extracted
1
  • As per the Tour there should be only one question asked per question.
    – PolyGeo
    Jul 20 '21 at 11:43
5

Do I need to set proj to NAD83 when defining the Point geometry too: point2 = ee.Geometry.Point([lon, lat],'EPSG:4269')

Yes, you must set this to correctly interpret the coordinates.

This is the only coordinate system that matters when retrieving point values; the crs in reduceRegion is irrelevant because the feature has no area — except for your next question.

How can I extract a 3x3 grid of pixels (9 sets of spectral values) centered at the point from Landsat instead of just a single (30mx30m) pixel at the lat-lon point?

One way to do this is to preprocess the image with neighborhoodToArray:

image.neighborhoodToArray(ee.Kernel.square(1))

Now each image pixel is a 2D array of the 3×3 group of values you want. However, the 3x3 area is computed over the output pixel grid, not the original image's, which means the CRS used by reduceRegion is now relevant. This CRS should not be the one of the point, but the one of the grid you want pixels from.

You can't just say “the pixels of the original image” because you're taking the median of a collection, and each image in the collection may have a different projection — even when they're from the same satellite, each orbit is slightly different. A useful approximation is to use the projection of some one of the images you're taking the median of. Or you can just choose a scale, but 3x3 patches are more likely to be sampled poorly if they don't match the image.

So, here's a working example (in JavaScript so it can be run in the EE Code Editor):

var point = ee.Geometry.Point([-121.771399, 40.330908],'EPSG:4269');

var landsatColl = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
    .filterBounds(point)
    .filterDate('2005-03-11', '2005-05-10');
var projection = landsatColl.first().projection();
var median = landsatColl.median();

var neighborhoods = median.neighborhoodToArray(ee.Kernel.square(1), 999);

var extracted = neighborhoods.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: point,
  crs: projection,
  scale: projection.nominalScale(),
});

Map.addLayer(median, {
  bands: ['B3', 'B2', 'B1'],
  min: 100,
  max: 2000,
  gamma: 1.4,
});
Map.addLayer(point);
print(extracted);

https://code.earthengine.google.com/c2819c5b28bc90b4186a8f29446767cd

This example using your point will show the point and pixels like this:

Pixel neighborhood

By using the Inspector and clicking on the pixel that the point occupies, we can see that the band values are B1: 404 B2: 550 B3: 594 B4: 2255 B5: 1573 B6: 2859 B7: 944. The printed results are:

B1: [[407,404,409],[426,404,429],[439,428,449]]
B2: [[570,570,590],[570,550,570],[590,571,571]]
B3: [[611,618,598],[578,594,579],[598,598,560]]
...

and the middle entry of these 3×3 nested lists (the downloaded form of arrays) give the same numbers.


There's one more thing you should do. You mention wanting several points; for this purpose, you should use reduceRegions instead of reduceRegion. The difference is that reduceRegions takes a FeatureCollection which can contain many points; this will be more efficient than calling reduceRegion for each point, because Earth Engine can query the image collection for images that intersect any of the points, rather than a separate query for each point.

4
  • Thanks a lot, @Kevin Reid. Super clear and helpful. Using your example, I can get 3x3 values. However, the spectral values I get by just intersecting the single point (and not the 3x3 kernel) do not match the values at the center of the 3x3 grid. In fact, they match the value right to the side of the center value (value at column 1 and row 2 of the 3x3 grid). Any idea what is happing here?
    – Mohamad
    Jul 19 '21 at 23:40
  • @Mohamad That's not what I see. I've added a description of the results I get to the answer — are you getting something different than that?
    – Kevin Reid
    Jul 20 '21 at 2:02
  • thanks @Kevin, you are absolutely right. My bad! I was using incorrect projection in reduceRegion. I can see the exact same results as you showed. On your last point about using reduceRegions, the problem is that I have a different date for each point which I use to get a time window (plus and minus one month) to compute median Landsat values among all images within that time window. Can I still use reduceRegions for a single query?
    – Mohamad
    Jul 20 '21 at 3:41
  • @Mohamad There might be a way to do it better, but not a straightforward one; go ahead and use reduceRegion.
    – Kevin Reid
    Jul 20 '21 at 4:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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