I am trying to plot a .tif raster on my map with basemap. With QGIS I see the raster layer as it is supposed to be: QGIS image
However, when then plotting with python basemap the colours are off, and the projection is somehow both rotated 180 degrees and mirrored, with a random blue line projected on the left side of the chart:
Some information about the .tif file (acquired with the rasterio and earthpy packages):
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7fcf9bdbef90> >
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 2760, 'height': 1350, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.01, 0.0, 3.7,
0.0, -0.01, 71.2)}
EPSG:4326
BoundingBox(left=3.7, bottom=57.7, right=31.3, top=71.2)
+proj=longlat +datum=WGS84 +no_defs
The original raster data was downloaded here (forest restoration potential) and cropped to the right latitude and longitude with QGIS.
What am I doing wrong?
import gdal
from numpy import linspace
from numpy import meshgrid
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
lowlong = 3.7 #lower left corner of longitude
lowlat = 57.7 #lower left corner of latitude
upplong = 31.3 #upper right corner of longitude
upplat = 71.2 #upper right corner of latitude
pathToRaster = r'~/Data/Shapefile/NorwayPotential.tif'
raster = gdal.Open(pathToRaster,1)
print(raster)
geo = raster.GetGeoTransform()
geo = raster.ReadAsArray()
mp = Basemap(projection='merc',
llcrnrlon=lowlong,
llcrnrlat=lowlat,
urcrnrlon=upplong,
urcrnrlat=upplat,
resolution='i')
mp.drawcoastlines()
mp.drawcountries()
x = linspace(0,mp.urcrnrx,geo.shape[1])
y = linspace(0,mp.urcrnry,geo.shape[0])
xx,yy = meshgrid(x,y)
mp.pcolormesh(xx,yy,geo)
plt.show()
After the line of code:
geo = raster.ReadAsArray()
you can flip the data array by
geo = geo[::-1,:]
and should get the correct result.