0

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:

Python basemap image

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()
Trenton McKinney
  • 56,955
  • 33
  • 144
  • 158
ElsRi
  • 37
  • 1
  • 5
  • Try using negative size pixel in the y direction – Ian Turton Sep 10 '20 at 19:02
  • How do I do that? – ElsRi Sep 11 '20 at 09:24
  • multiply the Y values you have by -1 ? – Ian Turton Sep 11 '20 at 09:30
  • If I do that the plot is entirely empty. Then I only see the basemap, not the data. I believe the y-values are coordinates, by making them negative the projection will be on the Southern half of the earth, in stead of in these countries – ElsRi Sep 11 '20 at 09:57
  • when you convert the Y values to pixels you should set the height of the pixel negative to allow for the difference in origins between images and geography – Ian Turton Sep 11 '20 at 10:48
  • Okay I found some additional explanation on the theory behind your comment. Sounds logical, but I am not very familiar with converting images to geography or gdal. How do I set the negative pixel height? Should I add something to the GetGeoTransform function? – ElsRi Sep 11 '20 at 11:15

1 Answers1

2

After the line of code:

geo = raster.ReadAsArray()

you can flip the data array by

geo =  geo[::-1,:]

and should get the correct result.

swatchai
  • 17,400
  • 3
  • 39
  • 58