1

I am currently working on a model that calculates the water balance of a basin "the VIC model". First we need to align (coregister) all of our input raster data sets to the grid template at the spatial resolution we need for all of the input datasets. The command we use is snap_raster.py which allows you to align the grid cells of two rasters and does resampling if needed:

Usage: python snap_raster.py <input raster> <output raster> <template raster> <subgrid boolean> <resampling method>

So this is the command:

Run commands: $ python snap_raster.py
../data/input/gis/Nyando_basin_SRTM_filled.tif
../data/input/gis/Nyando_basin_SRTMSnap.tif
../data/input/gis/Nyando_grid.tif True bilinear

But finally I don't get the output data but I get this message:

Traceback (most recent call last):
  File "snap_raster.py", line 143, in <module>
    main()
  File "snap_raster.py", line 137, in main
    snap_raster(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
  File "snap_raster.py", line 74, in snap_raster
    raise KeyError('Error setting output data type')
KeyError: 'Error setting output data type'

I'd like to know where the problem is and how to fix it.

Here is the snap_raster.py script:

# import dependencies
import sys
import numpy as np
from osgeo import gdal
from osgeo.gdalnumeric import *  
from osgeo.gdalconst import *

def snap_raster(inputRas,outputRas,templateRas,subGrid,resample):
    """
    FUNCTION: snap_raster
    ARGUMENTS: inputRas - input raster file to be snapped with TIF file extension
               outputRas - output snapped raster file with TIF file extension
               templateRas - input raster that the raster will be snapped to with TIF file extension
               subGrid - boolean value to set whether the output raster's 
                         resolution will be at the template raster resoltion or not
               resample - resampling method to be used when snapping
    KEYWORDS: n/a
    RETURNS: n/a
    NOTES: Returns no variables but writes an output file
    """

    # try to do the process
    try:

        # read source raster
        src = gdal.Open(inputRas,GA_ReadOnly)
        src_proj = src.GetProjection()
        srcXSize = src.RasterXSize
        srcYSize = src.RasterYSize

        # get no data value
        NoData = int(src.GetRasterBand(1).GetNoDataValue())

        # get source raster datatype
        dtype = gdal.GetDataTypeName(src.GetRasterBand(1).DataType)

        # get destination data type through a lookup dictionary
        datatypes = {'Byte':gdal.GDT_Byte,'Int16':gdal.GDT_Int16,
                    'UInt16':gdal.GDT_UInt16,'UInt32':gdal.GDT_UInt32,
                    'Int32':gdal.GDT_Int32,'Float32':gdal.GDT_Float32,
                    'Float64':gdal.GDT_Float64}

        src_dtype = datatypes[dtype]

        # try getting the resampling method        
        try:
            samps = {'nearest':GRA_NearestNeighbour,'bilinear':GRA_Bilinear,'cubic':GRA_Cubic,
                    'spline':GRA_CubicSpline, 'mean':GRA_Average, 'mode':GRA_Mode}  
            sampMethod = samps[resample]
        # if the resampling method does not exist, give an error
        except KeyError:
            raise KeyError('{0} is not a valid resampling method'.format(resample))

        # check to make sure that the output data type will make sense with the resampling method
        if (resample != 'nearest') | (resample != 'mode'):
            if 'Int' in dtype:
                src_dtype = datatypes['Float32']
            elif 'Byte' in dtype:
                src_dtype = datatypes['Float32']
            else:
                raise KeyError('Error setting output data type')

        # We want a section of source that matches this:
        match_ds = gdal.Open(templateRas,GA_ReadOnly)
        match_proj = match_ds.GetProjection()
        match_geotrans = match_ds.GetGeoTransform()
        matchXSize = match_ds.RasterXSize
        matchYSize = match_ds.RasterYSize

        # check if the subGrid argument is true or no
        if (subGrid == "True") | (subGrid == "true") | (subGrid == "TRUE"):
            # calculate the new raster resolution closest to source resolution
            xRatio = int(np.round(srcXSize / matchXSize))
            yRatio = int(np.round(srcYSize / matchYSize))

            wide = matchXSize * xRatio
            high = matchYSize * yRatio

            # create ouput geotransform with new raster resolution
            outGeom = [match_geotrans[0],(match_geotrans[1]/float(xRatio)) ,0,
                    match_geotrans[3],0,(match_geotrans[5]/float(yRatio))]

        else:
            # set output geotransform to the template raster
            wide = matchXSize
            high = matchYSize
            outGeom = match_geotrans

        # create output / destination raster
        dst = gdal.GetDriverByName('GTiff').Create(outputRas, wide, high, 1, src_dtype)
        dst.SetGeoTransform(outGeom)
        dst.SetProjection(match_proj)
        band = dst.GetRasterBand(1)
        band.SetNoDataValue(NoData)

        # Do the work
        gdal.ReprojectImage(src, dst, src_proj, match_proj, sampMethod)

        # flush the dataset
        dst = None 
        band = None

    # if not working, give error message
    except AttributeError:
        raise IOError('Raster file input error, check that all paths are correct')

    return

# main level program
def main():

    n_args = len(sys.argv)

    # Check user inputs
    if n_args != 6:
        print("Wrong user input")
        print("Script converts shapefile to raster grid")
        print("usage: python snap_raster.py <input raster> <output raster> <template raster> <sub grid boolean> <resampling method>")
        print("Exiting system...")
        sys.exit()

    else:
        # Pass command line arguments into function
        Snap_Raster(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])

    return

# Execute the main level program if run as standalone
if __name__ == "__main__":
    main()
1

You have a logic error here:

        else:
            raise KeyError('Error setting output data type')

Remove this.

Why?

If your datatype is int or byte and the resampling method is nearest or mode then an int/byte datatype is an appropriate output datatype. The code checks this, then casts the output datatype to Float32 for int/byte input with non-nearest/mode resampling. If the input is not int/byte it should just carry on as a float is an appropriate output data, but instead raises an error.

Also you have a logic error in:

if (resample != 'nearest') | (resample != 'mode'):

This will always return True as you need to test the resample isn't nearest AND isn't mode.. (also a minor point but can introduce subtle bugs, you are using bitwise not boolean operators. Use and and or instead of & and |).

E.g.

resample = 'nearest'
(resample != 'nearest') | (resample != 'mode')

True


resample = 'bilinear'
(resample != 'nearest') | (resample != 'mode')

True

It should be:

if (resample != 'nearest') and (resample != 'mode'):

E.g.

resample = 'bilinear'
(resample != 'nearest') and (resample != 'mode')

True

resample = 'nearest'
(resample != 'nearest') and (resample != 'mode')

False

Or even better:

resample not in ('nearest', 'mode')

So it should look like:

        # check to make sure that the output data type will make sense with the resampling method
        if resample not in ('nearest','mode'):
            if 'Int' in dtype:
                src_dtype = datatypes['Float32']
            elif 'Byte' in dtype:
                src_dtype = datatypes['Float32']

I would also recommend changing the bitwise operators here:

        if (subGrid == "True") | (subGrid == "true") | (subGrid == "TRUE"):

To or or subGrid in ('True', etc...) or even better:

        if subGrid.lower() == "True":

One last thing... Might be a typo, but your script shouldn't work at all. The 2nd last line of your main function calls Snap_Raster(sys.argv[1], etc...) but there is no Snap_Raster function defined, but there is a snap_raster. Python is case sensitive so Snap_Raster != snap_raster

  • Thank you very much for your answer it worked, and indeed I had to remove these logical errors. – Bouramtane Tarik Apr 9 '18 at 12:23

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.