I wrote function which allows me to plot raster and (optionally) plot vector on the same image. It works for small files, but in case of big rasters, I receive memmory error.
I've tried to add conditional statement to handle big rasters according to response in: Python/GDAL--Handling Big Rasters and avoid MemoryErrors?
However it doesn't work properly, I still receive an error message.
Could you help? Here is my Code:
def createMap(raster, vmax, vmin, output, shapefile=None, title=None):
###################################################################
# Creates image from raster and shapefile
# Based on: https://gist.github.com/jdherman/7434f7431d1cc0350dbe
######
# TODO: Consider rewriting to pyQGIS
# (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/composer.html)
#####
# Prerequisities:
# sudo apt-get install python-mpltoolkits.basemap
##################################################################
## Sample files for testing (comment in gedit: CTRL + M, uncomment: CTRL + SHIFT + M)
#import os
#from os.path import expanduser
#home = expanduser("~")
#SMOSfile = os.path.join(home,"Dropbox/Dane SMOS CATDS dla Wisły/DA_TC_MIR_CL_33/EXT-SM_RE02_MIR_CLF33A_20101231T000000_20120102T235959_272_001_7/ext-SM_RE02_MIR_CLF33A_20101231T000000_20120102T235959_272_001_7_1.DBL.nc")
#SMOSraster = 'NETCDF:"' + SMOSfile + '":Soil_Moisture'
#SentinelRaster = os.path.join(home,"Testy/calibrated_S1A_IW_GRDH_1SDV_20160512T161044_20160512T161.data/Sigma0_VH.img")
#vmin = 0
#vmax = 3000
#output = os.path.join(home,"testy.png")
#shapefile = os.path.join(home,"Dropbox/mapy/dorzecze_Wisły")
#createMap(SMOSraster, vmax, vmin, output, shapefile)
#createMap(SentinelRaster, vmax, vmin, output)
###################################################################
from osgeo import gdal, osr
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
# By default, osgeo.gdal returns None on error, and does not normally raise informative exceptions
gdal.UseExceptions()
gdata = gdal.Open(raster)
geo = gdata.GetGeoTransform()
xres = geo[1]
yres = geo[5]
# Get "natural" block size, and total raster XY size.
band = gdata.GetRasterBand(1)
block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
y_block_size = block_sizes[1]
xsize = band.XSize
ysize = band.YSize
print('x_block_size: {0}, y_block_size: {1}.'.format(x_block_size, y_block_size))
print('xsize: {0}, ysize: {1}.'.format(xsize, ysize))
if (xsize < 5000):
data = gdata.ReadAsArray()
else:
#########################################################
## TODO: for big rasters such as Sentinel-1:
## Solution adapted from https://gis.stackexchange.com/questions/211611/python-gdal-handling-big-rasters-and-avoid-memoryerrors
## It seems that I still receive Memory Error
y_block_size_NEW = int(round(y_block_size/200)) if y_block_size > 200 else y_block_size
x_block_size_NEW = int(round(x_block_size/200)) if x_block_size > 200 else x_block_size
# Create temporal raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(gdata.GetProjectionRef())
format = "GTiff"
driver = gdal.GetDriverByName( format )
# TODO: seems that I should create smaller temporal raster (?)
dst_ds = driver.Create("original_blocks.tif", xsize, ysize, 1, band.DataType )
dst_ds.SetGeoTransform(geo)
dst_ds.SetProjection(raster_srs.ExportToWkt())
blocks = 0
for y in xrange(0, ysize, y_block_size_NEW):
#print blocks
if y + y_block_size_NEW < ysize:
rows = y_block_size_NEW
else:
rows = ysize - y
for x in xrange(0, xsize, x_block_size_NEW):
if x + x_block_size_NEW < xsize:
cols = x_block_size_NEW
else:
cols = xsize - x
# Seems that some kind of average should be calculated here
array = band.ReadAsArray(x, y, cols, rows)
try:
array[array>0]=1
#print "we got them"
except:
print "could not find them"
dst_ds.GetRasterBand(1).WriteArray(array, x, y)
del array
blocks += 1
data = dst_ds.ReadAsArray()
# TODO: Remove temporal raster?
#########################################################
m = Basemap(llcrnrlon=17.00,llcrnrlat=48.75,urcrnrlon=25.25,urcrnrlat=54.50)
if shapefile is not None:
m.readshapefile(shapefile,'shp',drawbounds=True, color='0.3')
xmin = geo[0] + xres * 0.5
xmax = geo[0] + (xres * gdata.RasterXSize) - xres * 0.5
ymin = geo[3] + (yres * gdata.RasterYSize) + yres * 0.5
ymax = geo[3] - yres * 0.5
x,y = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
x,y = m(x,y)
cmap = plt.cm.gist_rainbow
cmap.set_under ('1.0')
cmap.set_bad('0.8')
im = m.pcolormesh(x,y, data.T, cmap=cmap, vmin=vmin, vmax=vmax)
cb = plt.colorbar( orientation='vertical', fraction=0.10, shrink=0.7)
if title is not None:
plt.title(title)
plt.savefig(output)