0

I have been using the Global Surface Water dataset for some analysis, but export the data to my local machine before using it. The dataset is unprojected on a WGS84 datum (EPSG 4326), but I want to analyze it in the native Landsat CRS, i.e. the respective UTM zone of each study site. When I use export, I'm simply entering the target CRS in the "Export" call. This is causing a couple of issues. When I compare the exported scene with the scene available in GSW, it looks like averaging between landsat 7-derived no data bands and existing water bodies (like channels) leads to spurious land pixels. In addition I notice significant warping of objects like lakes. An example comparison is below. The coordinates of this location are 69.601983 deg N, 161.742794 deg E.

This leads to the two questions I have:

1) Is specifying the CRS a reprojection from EPSG 4326 to EPSG X or redefining the CRS of the exported map (incorrectly)?

2) If this is a reprojection what is the reprojection method? The GSW data is categorical but takes on integers 0 to 2 (0 = No data, 1 = Land, 2 = Water), so a nearest neighbor resampling won't cause the above issue, but averaging or interpolation will.

The script I use to export over a given region is shown here. I should note that this region is on the edge of UTM 57 N and UTM 58 N so part of the warping may be because of that (but hopefully not). I have compared using 58 N in the export and get similar issues. I know that one solution to this may be to mark the 0 values with NA (or whatever JS equivalent) before exporting, but the significant warping is still bothersome.

UPDATE

I found that the reprojection in GEE appears significantly off from what I thought I was downloading. The script below was used to download data in EPSG:32657 (UTM57N). Alternatively, I downloaded data specifying no CRS and a scale of 30. This exports an image in EPSG:4326 with res = .000269 deg. I tested reprojecting to 15 m and 30 m resolution UTM57N scenes and found much better correspondence (visually). Four scenes are shown below which highlight the issue with the GEE export. I understand the issue is when I pass "30" as an argument to scale in export and the crs as the UTM zone, so to tie it back to Question 1: this is probably a reprojection but what is the new scale argument doing?

/// Visual parameter constants
// Black signifies no data 0, grey is land 1, white is water
var VIS_WATER_MASK = {
  min: 0.0,
  max: 2.0,
  palette: ['black', 'grey', 'white']
};



///RECURRENCE
//var date = "recurrence_07";

//var grab = ee.Image('JRC/GSW1_1/MonthlyRecurrence/monthly_recurrence_06');
print(grab)
/// Selection of individual months
//////////////////////////////////////////////// May
 var date = "2000_07";
 print(date);
 var filename = 'JRC/GSW1_1/MonthlyHistory/' + date;
 print(filename);

 var grab = ee.Image(filename);
 print(grab)
  var water1 = grab.select("water");
 Map.addLayer({
   eeObject: water1,  
   visParams: VIS_WATER_MASK,
   name: date
   });

Export.image.toDrive({
  image: grab,
  description: date,
  scale: 30,
  region: Kolyma,
  crs: 'EPSG:32657',
  folder: 'Kolyma_GSW_32657',
  maxPixels: 9e8
});

This is a GSW mask in EPSG 4326 (i.e. unprojected) with .000269 deg resolution, exported with scale = 30. This is about 15-m a pixel here (at 69 N). 4326

This is if you use the simple script I put above to export, setting CRS = 32657 (UTM 57N) and scale = 30. GEE57

This is when doing projection from 4326 to 32657 at 15 m spatial resolution (automatically set without specifying a resolution) using QGIS instead of GEE: GEE4326to57_15m

This is when doing projection from 4326 to 32657 at 30 m spatial resolution (prespecified target resolution) using QGIS instead of GEE: GEE4326to57_30m

2

Earth Engine generates reduced-resolution versions of each image when it is ingested. Part of the upload process involves specifying how the reduced-resolution pixels should be derived -- "mean" is the most common answer, but for categorical data like this "mode" is more appropriate.

I believe the problem here is that (a) this GSW data was inappropriately pyramided with "mean", and (b) because of the somewhat distorted pixels at the high latitudes Earth Engine is choosing to derive your 30m UTM57 pixels from 60m EPSG:4326 pixels.

You can avoid this issue by forcing Earth Engine to work with the full resolution data:

https://code.earthengine.google.com/6ff685f6d52e893701969e42acb60c36

(Note that if you zoom out here things will get very slow and/or run out of memory, since you'll be using a very large number of input pixels for each output pixel, but as long as you view or export at close to 30m you should be fine.)

| improve this answer | |
  • The results from QGIS are replicated when using the r stars package, which uses gdal I believe to do reprojection/warping. In your script, it looks like GEE's reprojection is exactly the issue. If you look at regions where black pixels (0, no data) and white (2, water) border, a fake set of grey (1, land) pixels is generated. Unless I'm sorely mistaken, nearest neighbor resampling shouldn't lead to that artifact. Setting scale = 15 in the projection removes that artifact. – Heymans Jun 9 at 23:54
  • Following up on previous comments R test, I exported a region with crs: 4326 and scale = 30. In R, using stars::st_warp(x, crs = 32657, cellsize = 30) where x is the scene loaded in as a stars object yields quite similar results to QGIS. – Heymans Jun 10 at 0:06
  • Sorry, those two comments were to the original answer. The updated answer looks great. – Heymans Jun 10 at 0:23

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.