2

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)
3
  • 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. Commented Nov 5, 2016 at 16:57
  • 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()?
    – matandked
    Commented Nov 6, 2016 at 10:22
  • 1
    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 Commented Nov 6, 2016 at 14:02

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.