This is a rather conceptual question that i haven't been able to address yet. Let's say I have a raster file projected in an equal area projection (say Mollweide), with a cell resolution of 100km x 100km. Since the projection preserves the area along the parallels, all cells should have the same area.
Now, as we know from the basic principles of cartography, no projection can preserve both area and shape of features when converting a sphere to a flat surface. This can be easily illustrated, if we draw a line of 100 km projected length at the equator on a map in equal area projection (but also on a map in orthographic projection), and a line of the same projected length towards the poles, these two lines will appear to have very different lengths in terms of screen pixels (check it out in google earth!).
Given this knowledge, if we try to visualise the above mentioned raster either in Python or in R using the appropriate packages (see scripts below), the cells have the same screen pixel dimensions regardless of the latitude. However, I would expect that, since each cell has the same dimensions (100km x 100km), cells towards the poles should be thinner and rather elongated (see also the attached snapshot from a recent scientific publication) in order to preserve area. So my question is, what is the appropriate way of visualising raster cells in equal area projections?
Code in Python:
import matplotlib.pyplot as plt from osgeo import gdal import cartopy import numpy as np fig = plt.figure(figsize=(15, 8)) proj = cartopy.crs.Mollweide() ax = plt.axes(projection = proj) raster = gdal.Open('test.tif') # a raster at with cell resolution 100x100km band = raster.GetRasterBand(1) nodata = band.GetNoDataValue() array = band.ReadAsArray() array[array == nodata] = np.nan img = ax.imshow(array, origin='upper', cmap='jet') ax.get_xaxis().set_visible(True) ax.get_yaxis().set_visible(True) plt.show()
Code in R:
library(raster) ras = raster('test.tif') plot(ras)
Example figure with stretched cells towards the north pole (source):