1

I would like to divide the moderate resolution satellite imagery into smaller pixels. To clarify this, I took a Landsat 30x30m image. For further analysis, I first want to split all pixels into 3x3m pixels, exactly corresponding with the original pixels. The method ee.Image.reproject() seems to come closest, although the pixels are reprojected into a slight different projection and therefore the pixels do not match with the original Landsat pixels. How can I Create an image where all new smaller pixels exactly match the original image?

Here a piece of simplified code which shows the mismatch (zoomed in to the corner of the image, where the mismatch is most clearly seen):

The mismatch removes when '.clip' is done at the end, where the image is added to the map

// Extent for the images to clip
var Haurvig = ee.Geometry.Rectangle([8.1247, 55.927, 8.1625, 55.907]);

// Load an LS8 image with least clouds
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                     .filterDate('2018-01-01','2018-12-01')
                     .filterBounds(Haurvig)
                     .sort('CLOUD_COVER').first();
// rename the band names and clip the image
var image =  image.select(['B2','B3','B4','B5','B6','B7','B10','pixel_qa'])
          .rename(['blue','green','red','nir','swir1','swir2','tir','qa'])
          .clip(Haurvig);    

// Reproject the image to 3x3 m
var reprojected = image.reproject(image.projection().crs(), null, 3);

// add to the map
var trueColor = {bands: ['red','green','blue'], min: 300, max: 3800};
Map.addLayer(reprojected, trueColor, 'Reprojected image');
Map.addLayer(image, trueColor, 'Original image');
Map.centerObject(ee.Geometry.Point([8.162408851233522,55.907053640896]), 20);
  • Although what you are asking for is possible, it is often not the most efficient way to use Earth Engine. Stepping back a level, what are you planning on doing with the image with 3x3m pixels (i.e. what is your overall goal)? This may help others suggest a good solution. – Tyler Erickson Oct 20 '18 at 22:04
  • I didn't want to add that in my question, as that could start a discussion. However, what I would like to approach is the method described in this article (doi.org/10.1080/01431160500213292; super resolution mapping). Using first Linear spectral unmixing on a per-pixel basis, I would like to classify water, sand and vegetation on a sub-pixel scale, maintaining the original information on the pixels itself. The ultimate goal is to get a polygons/lines which shed purely water from the beach and marks the abrupt change from sand overblown dunes to vegetated dunes. – Kuik Oct 21 '18 at 7:46
1

I'd like to start off this answer by qualifying it to say that the approach described is useful for understanding how projections work, but it is rarely the best approach for large-scale image analysis workflows. The use of ee.Image.reproject() short-circuit's Earth Engine's default behavior of reprojecting as needed and using image pyramids to efficiently query data. Make sure to read the Earth Engine docs on Projections and Resampling to understand how ee.Image.reproject() works before using it.

That being said, here is a way to create a small sample image with pixels that are aligned with the original image.

//  Define a small Region of Interest (ROI) used to clip results.
var roi = ee.Geometry.Point(8.3431, 55.8825).buffer(60);

// Load an LS8 image with least clouds.
var image_1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                     .filterDate('2018-01-01','2018-12-01')
                     .filterBounds(roi)
                     .sort('CLOUD_COVER').first();

Get the projection information from the image, and scale it to create a new projection that has 10 times as many pixels in each dimension.

// Resample and reproject the image to the new projection.
var projection_1 = image_1.projection();
var projection_2 = projection_1.scale(1/10, 1/10);
var image_2 = image_1.resample('bicubic').reproject(projection_2);

Note that the .resample() method is used to make it easier to compare the original and reprojected pixels. Remove it if you want to use the default 'nearest neighbor' resampling method.

// Print out the projection information, for easy inspection.
print('projection_1', projection_1);
print('projection_2', projection_2);

Compare the two projections using a split map view.

// Display the various images as map layers on a split map. 
var TRUE_COLOR_VIS = {bands: ['B4', 'B3', 'B2'], min: 300, max: 600};
var leftMap = ui.Map();
var rightMap = ui.Map();
leftMap.addLayer(image_1.clip(roi), TRUE_COLOR_VIS, 'image_1 (original image)');
rightMap.addLayer(image_2.clip(roi), TRUE_COLOR_VIS, 'image_2 (reprojected to 3m postings)');

// Create a SplitPanel to hold the adjacent, linked maps.
var splitPanel = ui.SplitPanel({
  firstPanel: leftMap,
  secondPanel: rightMap,
  wipe: true
});
var linker = ui.Map.Linker([leftMap, rightMap]);

// Set the SplitPanel as the only thing in the UI root.
ui.root.widgets().reset([splitPanel]);
leftMap.centerObject(roi);

Here is an interactive app for comparing the projections: https://tylere.users.earthengine.app/view/gis-stackexchange-299589

enter image description here

Source code is here: https://code.earthengine.google.com/e7f2fa092999faec3209d07dfd41d5dd#

Note in the original code, the "mismatch" of pixels along the edge was caused by the relative placement of the .clip() method, which occurred before the .reproject() method. In general, it is a good practice to only use .clip() at the end of your analysis (i.e. when displaying results on the map).

  • Thanks a lot. The only thing missing in my code was thus that I needed to clip in the end of my analysis. – Kuik Oct 23 '18 at 10:52
1

You could reduce pixel size with gdalwarp command outside of google-earth-engine and, afterward, upload there modified raster. Following qgis python code (pyqgis) works well.

import os

raster = iface.activeLayer()
provider = raster.dataProvider()

crs = raster.crs().authid()

input_raster_path = provider.dataSourceUri()

output_raster_path = '/home/zeito/pyqgis_data/NDVI/b3.ND_clip_3_3_res.tif'

root, filename = os.path.split(output_raster_path)

cmd = "gdalwarp -overwrite -t_srs " + crs + " -tr " + str(xsize) + " " + str(ysize) + " -s_srs " + crs + " " + \
      input_raster_path + " " + output_raster_path

os.system(cmd)

changed = QgsRasterLayer(output_raster_path,
                         filename)

QgsMapLayerRegistry.instance().addMapLayer(changed)

I tried it out with raster (30x30 m) of following image:

enter image description here

Resulting raster has expected resolution (3x3 m) and it is perfectly aligned with original raster.

  • Although I am sure it would work fine, the purpose of my research is to develop a systematic way to work this out inside the GEE. Unfortunately, this is not the solution for me. If there isn't any proper way to do it inside the GEE, I will need to find another approach to coop with the sub-pixel classification of the LS imagery. – Kuik Oct 22 '18 at 9:03

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.