8

I am trying to plot raster map using Pyhton. Image is not being saved, though it appears. Moreover, pixels are not appearing on the exact place. I am getting image as shown below (which is WRONG, however, lat-lon labeling is correct.)

enter image description here

Below is the correct orientation, which I have exported from ArcGIS. enter image description here

There is some mistake in placing the pixel at right location, but I ma not getting it. How to correct it?

Here is my script.

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from numpy import linspace
from numpy import meshgrid
import numpy as 

#################################################################################

ds = gdal.Open('Path\\To\\Raster.tif')
data = ds.ReadAsArray()

gt = ds.GetGeoTransform()   
proj = ds.GetProjection()

#################################################################################

xres = gt[1]
yres = gt[5]

xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5

x_center=(xmin+xmax)/2
y_center=(ymin+ymax)/2    
#################################################################################

fig = plt.figure(figsize=(18, 10))
ax = fig.add_subplot(111, axisbg='w', frame_on=True)

m = Basemap(llcrnrlon=xmin,llcrnrlat=ymin,urcrnrlon=xmax,urcrnrlat=ymax,
             projection='tmerc', lat_0 = y_center, lon_0 = x_center)

parallels = np.arange(15.,30.,0.25)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=12, linewidth=0.4) 

meridians = np.arange(70.,90.,0.25)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=12, linewidth=0.4)

x = linspace(0, m.urcrnrx, data.shape[1])
y = linspace(0, m.urcrnry, data.shape[0])

xx, yy = meshgrid(x, y)

m.pcolormesh(xx, yy, data, cmap=plt.cm.jet)

plt.show()
plt.savefig('Path\\To\\Save_Image.png', bbox_inches='tight', pad_inches=.2, dpi=600)    
  • Have you tried imshow instead of pcolormesh? It takes a different corner as the image origin by default. – countryman Mar 13 '16 at 18:00
  • When I replace pcolormesh with imshow it says NameError: name 'imshow' is not defined How to use `imshow' here? @countryman – Mario Mar 13 '16 at 18:10
  • Here are some really good examples: matplotlib.org/basemap/users/examples.html – countryman Mar 13 '16 at 19:05
  • 1
    it is flipped up down have you flipped the origin of the raster which can be either bottom left or top left in some applications – user681 Mar 13 '16 at 19:40
  • @DanPatterson I have no idea, please have a look into my script to find where I am messing it. – Mario Mar 14 '16 at 4:40
6

I solved this problem.

I was reading the raster file row-wise from top to bottom, and plotting it row-wise from bottom to top.

Since, I can't do anything with plotting I flipped the array, and it worked.

Below is the modification which I made,

ds = gdal.Open('Path\\To\\Raster.tif')
data = ds.ReadAsArray()
data = np.flipud(data)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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