Create map (image) from raster and vector in Python - problem with big files

Question

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)

• user6072577: It seems like the else block is splitting the data into separate chunks before putting it back together again, so that nothing has actually changed by the time "dst_ds.ReadAsArray()" is called. Could you rewrite it so that multiple datasets are created and then plot these individually? Also making the blocksize 1/200 of the original is probably excessive unless the dataset is truly massive, maybe it would be better to split the image in 2*2 or 3*3 blocks.
• matandked: Thank you for suggestion. Will plotting of those pieces separately not overwrite previously plotted piece? I'm not sure if I also need such level of details - will it not be easier to generalize/simplify raster, so that it will be easier to ReadAsArray()?
• user6072577: If you don't need the full resolution then the easiest solution will be to resample the image to a lower resolution. You can do this with gdal warp from the command line: gdalwarp -of gtiff -tr 100 100 input_name.tif output_name.tif