Trazar el ráster GDAL usando matplotlib Basemap

Me gustaría trazar un raster tiff ( descargar -723Kb) usando el mapa base matplotlib. Las coordenadas de proyección de mi raster es en metros:

In [2]: path = r'albers_5km.tif' raster = gdal.Open(path, gdal.GA_ReadOnly) array = raster.GetRasterBand(20).ReadAsArray() print ('Raster Projection:\n', raster.GetProjection()) print ('Raster GeoTransform:\n', raster.GetGeoTransform()) Out [2]: Raster Projection: PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",15],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",30],PARAMETER["longitude_of_center",95],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]] Raster GeoTransform: (190425.8243, 5000.0, 0.0, 1500257.0112, 0.0, -5000.0) 

Si trato de trazar esto usando una proyección de Robin usando contourf con latlon=False se supone que las coordenadas de proyección de mapas son coordenadas de proyección de mapa (ver documentos , creo que eso es lo que tengo).

Pero si miro la ttwig, me doy cuenta de que está situada abajo a la izquierda muy pequeña:

abajo a la izquierda

Utilizando este código:

 In [3]: xy = raster.GetGeoTransform() x = raster.RasterXSize y = raster.RasterYSize lon_start = xy[0] lon_stop = x*xy[1]+xy[0] lon_step = xy[1] lat_start = xy[3] lat_stop = y*xy[5]+xy[3] lat_step = xy[5] fig = plt.figure(figsize=(16,10)) map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0) lons = np.arange(lon_start, lon_stop, lon_step) lats = np.arange(lat_start, lat_stop, lat_step) xx, yy = np.meshgrid(lons,lats) levels = [array.min(),-0.128305,array.max()] map.contourf(xx, yy,array, levels, cmap=cm.RdBu_r, latlon=False) map.colorbar(cntr,location='right',pad='10%') map.drawcoastlines(linewidth=.5) map.drawcountries(color='red') 

Eventualmente no quiero tener una vista del mundo sino una vista detallada. Pero esto me da un nivel de zoom donde se dibujan las líneas costeras y los países, pero los datos se colocan de nuevo en la esquina inferior izquierda, pero no tan pequeños como en el tiempo anterior:

de nuevo abajo a la izquierda

Usando el siguiente código:

 In [4]: extent = [ xy[0],xy[0]+x*xy[1], xy[3],xy[3]+y*xy[5]] width_x = (extent[1]-extent[0])*10 height_y = (extent[2]-extent[3])*10 fig = plt.figure(figsize=(16,10)) map = Basemap(projection='stere', resolution='c', width = width_x , height = height_y, lat_0=40.2,lon_0=99.6,) xx, yy = np.meshgrid(lons,lats) levels = [array.min(),-0.128305,array.max()] map.contourf(xx, yy, array, levels, cmap=cm.RdBu_r, latlon=False) map.drawcoastlines(linewidth=.5) map.drawcountries(color='red') 

Puede usar el siguiente código para convertir las coordenadas, toma automáticamente la proyección de su ráster como origen y la proyección de su objeto de mapa base como el sistema de coordenadas objective.

Importaciones

 from mpl_toolkits.basemap import Basemap import osr, gdal import matplotlib.pyplot as plt import numpy as np 

Conversión de coordenadas

 def convertXY(xy_source, inproj, outproj): # function to convert coordinates shape = xy_source[0,:,:].shape size = xy_source[0,:,:].size # the ct object takes and returns pairs of x,y, not 2d grids # so the the grid needs to be reshaped (flattened) and back. ct = osr.CoordinateTransformation(inproj, outproj) xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T)) xx = xy_target[:,0].reshape(shape) yy = xy_target[:,1].reshape(shape) return xx, yy 

Lectura y procesamiento de los datos.

 # Read the data and metadata ds = gdal.Open(r'albers_5km.tif') data = ds.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 # create a grid of xy coordinates in the original projection xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres] 

Trazar

 # Create the figure and basemap object fig = plt.figure(figsize=(12, 6)) m = Basemap(projection='robin', lon_0=0, resolution='c') # Create the projection objects for the convertion # original (Albers) inproj = osr.SpatialReference() inproj.ImportFromWkt(proj) # Get the target projection from the basemap object outproj = osr.SpatialReference() outproj.ImportFromProj4(m.proj4string) # Convert from source projection to basemap projection xx, yy = convertXY(xy_source, inproj, outproj) # plot the data (first layer) im1 = m.pcolormesh(xx, yy, data[0,:,:].T, cmap=plt.cm.jet) # annotate m.drawcountries() m.drawcoastlines(linewidth=.5) plt.savefig('world.png',dpi=75) 

introduzca la descripción de la imagen aquí

Si necesita que la ubicación de los píxeles sea 100% correcta, es posible que desee verificar la creación de matrices de coordenadas con mucho cuidado (porque no lo hice en absoluto). Esperamos que este ejemplo lo ponga en el camino correcto.