World Settlement Footprint - Binary variable changing to continuous variable after coordinate referencing system change

Question

I am working with World Settlement Footprint data from 2019, which uses Sentinel data to create a binary mask of settlement area (=255 if settlement, =0 otherwise). I encountered a strange result when switching projections in R, and I was wondering if anyone can explain what is causing the issue and if there is any way around it (either directly in R, or possibly just exporting from Earth Engine as my desired projection).

To get the data in my desired format, I first downloaded a geotiff of one tile from this website https://download.geoservice.dlr.de/WSF2019/ then uploaded it to Earth Engine.

Then, I exported as a geotif using the default CRS (EPSG:3857). My script can be accessed here:https://code.earthengine.google.com/bf17e647c502d101c19303f3bfe588e6

Uploading my data into R and mapping it I get unique values of 0 and 255, which is to be expected.

However, when I change the CRS, the unique values are all over the place, and sometimes even go negative (values seem to be roughly in the range of [-40,300]). Visually the settlement map looks similar (i.e. more densely settled areas are still denser), but I'm not sure why changing the CRS is assigning values other than 0 and 255 to different pixels. Seemingly my code indicates that there are now over 200 thousand unique values (while previously there were only 2)

#World settlement footprint. Binary indicator of whether an area is a settlement or not
wsf <- raster("Data/10m_WSF_2019_cb.tif")

# UTM zone for Cox's Bazar ( zone 46N). This will be used later to accurately measure euclidian distances over the area.
crs_utm = "+proj=utm +zone=46 +datum=WGS84 +units=m"


# Project the raster to UTM
wsf<- projectRaster(wsf, crs = crs_utm)
plot(wsf)
   
#print unique values
unique(values(wsf))
# [1]   [1]            NA   0.000000000  50.613453125 181.947251063 112.249735884 -67.419255612 -33.270774722
# [8] -69.337822817 -33.962599669 172.938216558 255.000000000 255.000000000   1.029044815 -37.617426678........ [ reached getOption("max.print") -- omitted 210415 entries ]

Answer 1

reprojecting a raster uses by default bilinear interpolation:

 projectRaster(from, to, res, crs, method="bilinear", 
              alignOnly=FALSE, over=FALSE, filename="", ...) 

because each cell in the output doesn't overlap exactly the inputs cells (because of the distortions in projections) this interpolation causes cells in categorical rasters to have non-categorical values.

So what method to use?

 method: method used to compute values for the new RasterLayer. Either
          'ngb' (nearest neighbor), which is useful for categorical
          variables, or 'bilinear' (bilinear interpolation; the default
          value), which is appropriate for continuous variables.

"ngb", for nearest neighbour, most likely.

Also, obligatory "try using the terra package instead of raster if you can" message. Functions have generally different names but raster is basically frozen now and terra is the future.

Answer 2

You can do

library(terra)
url <- "https://download.geoservice.dlr.de/WSF2019/files//WSF2019_v1_92_20.tif"
x <- rast(url)
to_crs <- "+proj=utm +zone=46 +datum=WGS84 +units=m"
y <- project(x, to_crs, method="near", threads=TRUE, 
             filename="WSF2019_v1_92_20_utm.tif", datatype="INT1U")

I am not sure if this makes sense to do if your purpose is "to accurately measure Euclidian distances over the area.". Transforming lon/lat data to a planar coordinate reference system makes distance computation simpler (and perhaps that is your purpose), but generally, only less accurate measurements can be made because of the distortion introduced.