2

I confused the setting of Basemap about Alers projection. I want to plot a province area of Nation to pcolormesh. the value of lat_0 and lon_0 is Nation map's center longitude and latitude or this province map's center and origin lat?

driver=gdal.GetDriverByName('GTiff')
driver.Register()
file0=r'I:\extract\temp\m2000041.tif'
ds=gdal.Open(file0,gdal.GA_ReadOnly)
band=ds.GetRasterBand(1)
arr=band.ReadAsArray()
gt=ds.GetGeoTransform()
proj=ds.GetProjection
xres = gt[1]
yres = gt[5]

# get the edge coordinates and add half the resolution 
# to go to center coordinates
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

ds = None
width = xmax - xmin + xres
height = ymax - ymin - yres


lllon, lllat, urlon, urlat = gt[0], ymin - 0.5 * yres, xmax + 0.5 * xres,  ymax - 0.5 * yres
xl, yl = np.arange(lllon,urlon,xres), np.arange(urlat,lllat,yres)

false_easting = 4000000                 #False easting
xm,ym = np.meshgrid(xl-false_easting,yl)
arr2 = arr
arr2[arr2==-32767]=0                 #exclude nodata value

fig = plt.figure()                 
plt.title("The Plot")
m = Basemap(width = width,height = height,projection = 'aea',lat_1 = 25,lat_2 = 47,lon_0 = 105,lat_0 = 0,resolution = 'c')  # set the Albers paraller lat1,lat2, central lon and lat origin of Nation map.
lonl,latl = m(xm,ym)

m.pcolormesh(lonl,latl, arr2)
plt.show()

This code output a blank plot. I doubt the coordiantion meshgrid is wrong. How to plot the Albers map? Would you give a example?

2

Given a projected raster written to a tif file, I use the following code snippet to plot it with Basemap. It uses rasterio to read the data and pyproj for coordinates transformations.

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from pyproj import Proj

file = '/path/to/file.tif'

# Read the file data and metadata
with rasterio.open(file) as src:
    meta = src.meta
    data = src.read(1, masked=True)

# Extract shape/extent of the raster in its CRS
map_width = meta['width'] * meta['affine'][0]
map_height = meta['height'] * meta['affine'][0]
xmin = meta['affine'][2]
xmax = xmin + map_width
ymax = meta['affine'][5]
ymin = ymax - map_height
llproj = (xmin, ymin)
urproj = (xmax, ymax)
extent = [xmin, xmax, ymin, ymax] # [left, right, bottom, top]

# Instantiate projection class and compute longlat coordinates of
# the raster's ll and ur corners
crs = meta['crs']
p = Proj(**crs)
llll = p(*llproj, inverse=True)
urll = p(*urproj, inverse=True)

# Instantiate Basemap class
m = Basemap(llcrnrlon=llll[0], llcrnrlat=llll[1], urcrnrlon=urll[0], urcrnrlat=urll[1],
            projection=crs['proj'],
            resolution='l', lat_0=crs['lat_0'], lon_0=crs['lon_0'])
            # There might be other parameters to set depending on your CRS


# draw coastlines.
m.drawcoastlines()
# PLot raster
m.imshow(data, origin='upper', extent = extent)
m.colorbar()
# draw parallels and meridians.
parallels = np.arange(0.,81,5.)
meridians = np.arange(10.,351.,5.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[True,False,False,False])
m.drawmeridians(meridians,labels=[False,False,False,True])
plt.show()
plt.savefig('figure_name.png',dpi=100, transparent=True)
  • This is great. Very useful snippet! – Tactopoda Feb 11 '18 at 0:19

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.