1

I have a raster projected in EPSG:21781 (available here) that I'd like to display with mapview. Raster values correspond to the grid id, and need to be preserved. If I understood correctly, mapview is calling leaflet::addRasterImage() to do so, and reproject the raster into EPSG:3857 (pseudo-mercator) through raster::projectRaster(). The first problem is that the reprojection is done through a 'bilinear' method (by default) whereas I need the nearest-neighbor method. To overcome this problem, I have reprojected the raster manually using the following command:

grid <- raster::raster("~/grid100.tif")

grid2 <- projectRaster(grid,crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs",method="ngb") # corresponding to epsg 3857

summary(values(grid))

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1035406 2070572 2074309 3111142 4157490 3554337

summary(values(grid2))

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1056410 2097896 2093243 3133844 4157490 3917853

table(is.na(match(values(grid),values(grid2))))

FALSE TRUE
6206341 1488539

Conclusion: 1'488'539 cells are dropped after reprojection... Any idea of what is going on here? Thanks

Improve this question
6
  • 3
    Whenever you reproject a raster from one coordinate system to another, you are regridding the data. Between some projections, that regridding will lead to changes in number of pixels in the resulting raster. This is the general reason for avoiding reprojecting the data as much as possible.
    – Mikkel Lydholm Rasmussen
    Commented Feb 28, 2020 at 14:20
  • Thanks for replying. The issue here is that mapview or leaflet "automatically" reproject rasters into EPSG:3857 (for compatibility with basemap providers). Does it mean that if one want to keep the resolution of a raster, one can't use those packages?
    – Ervan
    Commented Feb 28, 2020 at 15:13
  • 1
    If you want to keep the exact resolution, yes. But the reprojection is not that bad. Right now, you are directly comparing the values, which are likely changed slightly by the reprojection. If you instead look at the total number of pixels, you'll notice a much smaller difference between the layers.
    – Mikkel Lydholm Rasmussen
    Commented Feb 28, 2020 at 15:30
  • 1
    +1 to the above. It's not that any cells have been dropped, they just have a different value assigned to them. The number of cells should still be the same: ncell(grid) and ncell(grid2) should still be identical.
    – Where's my towel
    Commented Mar 2, 2020 at 13:59
  • 1
    @Where'smytowel: ncell(grid) = 7694880 while ncell(grid2)= 7751910
    – Ervan
    Commented Mar 3, 2020 at 9:32

1 Answer 1

Reset to default
2

I think I found the problem: for some reason, the projectRaster command introduces a "frame" of additional, empty (NA) cells around the original data.

As discussed in the comments above, there are no cells being dropped, there are cells being added: ncell(grid) = 7694880 while ncell(grid2)= 7751910

The difference in the number of cells between the two grids is exactly the difference in the number of NA values:

ncell(grid) - ncell(grid2)
sum(is.na(grid[])) - sum(is.na(grid2[]))

Result is in both cases -57030, which is a hint that all new cells are likely to be empty/NA.

Investigating further, there are 10 rows and 10 columns more than before:

nrow(grid) - nrow(grid2) # = -10
ncol(grid) - ncol(grid2) # = -10

If you export grid2, throw both of them into a GIS and color the NA's differently, you can easily see that grid2 has some sort of 5 cells wide border. Alternatively, here's a quick and dirty (and ugly) illustrative plot in R:

# define a area at the top edge to zoom in to:
ext <- new("Extent", xmin = 682940.562171853, xmax = 692881.65578179, 
    ymin = 290620.553207836, ymax = 301271.724932769)

# see where it is:
plot(grid2)
plot(ext, add=TRUE, col="red")

# the actual plot:
plot(ext, type="n") # empty plot zoomed in to an area
plot(grid2, add=T, colNA="blue", maxpixels=ncell(grid2)) # add grid2, paint NA's blue, use "maxpixels" to force full resolution - might take a while to plot
plot(grid, add=T, colNA="grey", maxpixels=ncell(grid)) # add grid with grey NA's

The blue area contains the new NA's. I don't know why this happens, but it is easily solved with a simple crop:

grid_fixed <- crop(grid2, grid)
Improve this answer
4
  • 1
    Thanks for trying to push this forward. My problem isn't the creation of NA values (this can easily be solved), but rather that 'grid.id's that are present in grid aren't in grid2... table(is.na(match(values(grid),values(grid2)))) shows that 1'488'539 grid id 's don't match between the 2 rasters.
    – Ervan
    Commented Mar 7, 2020 at 13:47
  • 1
    Because the raster is interpolated when it is being projected, you can't really avoid that some values change. It doesn't matter if you use bi-linear or nearest neighbor interpolation, the values will change.
    – JonasV
    Commented Mar 7, 2020 at 16:12
  • Once you have the same number of rows and columns in both grids, can't you just copy the values from the old to the new? values(grid2) <- values(grid) should work afaict.
    – Where's my towel
    Commented Mar 9, 2020 at 5:54
  • Thanks for your inputs. crop() doesn't work here as rasters don't overlap (getting an error message). So cropping and reassigning values isn't possible. Not much one can do I am afraid.
    – Ervan
    Commented Mar 9, 2020 at 9:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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