9

I want to average a 30m raster at the resolution of a 1000m one. Specifically, I want to computer, for each Daymet raster cell, the percentage of cultivated land from the 30m CDL.

How can I do that in EarthEngine? My question is very similar to one asked on the EE list. The author was exploring 3 solutions:

  1. A bilinear resample
  2. Reproject
  3. ReduceRegions, using a coarser raster that has been converted to features.

Someone gave an answer using reduceNeighborhood() to use a mean reducer, then reproject(), but it does not give quite the right result (see below, not right scale). How can I adjust that solution (I'm giving wrong scale parameters presumably?), or is there another solution that would be better?

Code

// DATA Import
var CDL = ee.ImageCollection("USDA/NASS/CDL")
var daymet = ee.ImageCollection("NASA/ORNL/DAYMET_V3"). 
  filter(ee.Filter.calendarRange(2015, 2015, "year"))

var aoi = ee.Geometry.Polygon([[[-100.501, 42.9819],[-100.501, 42.4045], [-98.39, 42.40457], [-98.39, 42.98192]]])    

var CDL_2015 = ee.Image("USDA/NASS/CDL/2015").select("cultivated").clip(aoi)
var DYM_2015 = ee.Image(daymet.first()).select("tmin").clip(aoi)


// Operation
var image_frac=CDL_2015.eq(2).reduceNeighborhood({
  reducer: ee.Reducer.mean(),
  kernel: ee.Kernel.square(250,"meters"),
}).reproject(DYM_2015.projection().atScale(1000)).rename("cultivated")


//Visualization
Map.centerObject(aoi, 12)

Map.addLayer(DYM_2015.randomVisualizer(), {},  'DAYMET tmin')
Map.addLayer(CDL_2015.select("cultivated").eq(2), {min:0, max:1, opacity: 0.8, palette: ["beaed4","7fc97f"]}, "CDL coverage original")
Map.addLayer(image_frac.select("cultivated"), {min:0, max:1, opacity: 0.4}, "CDL coverage")
Matifou
  • 2,011
  • 1
  • 21
  • 43
  • The scale and the kernel radius seem to be required to be the same. Still, it does produce some erroneously zeros at location where there is little cover of the CDL2015 image. – Kuik Apr 17 '19 at 08:29

1 Answers1

3

Based on the sample (reproduced below) at https://developers.google.com/earth-engine/resample - I think that the following code should work where your existing reduceNeighborhood code is:

var image_frac=CDL_2015.eq(2).reduceResolution({
  reducer: ee.Reducer.mean(),
  maxPixels: 4096,  // something large enough to not cause it to error out
}).reproject(DYM_2015.projection()).rename("cultivated")

Here's their code and some notes from them:

// Load a MODIS EVI image.
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A1').first())
    .select('EVI');

// Display the EVI image near La Honda, California.
Map.setCenter(-122.3616, 37.5331, 12);
Map.addLayer(modis, {min: 2000, max: 5000}, 'MODIS EVI');


// Get information about the MODIS projection.
var modisProjection = modis.projection();
print('MODIS projection:', modisProjection);


// Load and display forest cover data at 30 meters resolution.
var forest = ee.Image('UMD/hansen/global_forest_change_2015')
    .select('treecover2000');
Map.addLayer(forest, {max: 80}, 'forest cover 30 m');


// Get the forest cover data at MODIS scale and projection.
var forestMean = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });


// Display the aggregated, reprojected forest cover data.
Map.addLayer(forestMean, {max: 80}, 'forest cover at MODIS scale');

They also note:

In this example, note that the output projection is explicitly set with reproject(). During the reprojection to the MODIS sinusoidal projection, rather than resampling, the smaller pixels are aggregated with the specified reducer (ee.Reducer.mean() in the example). This sequence of operations is illustrated in Figure 3. Although this example uses reproject() to help visualize the effect of reduceResolution(), most scripts don't need to explicitly reproject; see the warning here.

Note that a second reprojection occurs (implicitly) to display the data on the Code Editor map. Visually inspect the results and observe the correspondence between the pixels from the MODIS layer and the forest cover data reprojected to MODIS scale and projection. In general, you should rarely need to explicitly reproject() in Earth Engine.

Community
  • 1
nicksan
  • 1,716
  • 1
  • 12
  • 22
  • 1
    Also, make sure to specify the scale when reprojecting, as mentioned in the answer here: https://gis.stackexchange.com/a/356795/33177 – shahryar Apr 01 '20 at 19:21