geolib.py 文件源码

def gdaldem_wrapper(fn, product='hs', returnma=True, verbose=True):
    """Wrapper for gdaldem functions

    Note: gdaldem is directly avaialable through API as of GDAL v2.1

    https://trac.osgeo.org/gdal/wiki/rfc59.1_utilities_as_a_library

    This function is no longer necessry, and will eventually be removed.
    """
    #These gdaldem functions should be able to ingest masked array
    #Just write out temporary file, or maybe mem vrt?
    valid_opt = ['hillshade', 'hs', 'slope', 'aspect', 'color-relief', 'TRI', 'TPI', 'roughness']
    try:
        open(fn)
    except IOError:
        print("Unable to open %s" %fn)

    if product not in valid_opt: 
        print("Invalid gdaldem option specified")

    import subprocess
    from pygeotools.lib import iolib
    bma = None
    opts = []
    if product == 'hs' or product == 'hillshade':
        product = 'hillshade'
        #opts = ['-compute_edges',]
        out_fn = os.path.splitext(fn)[0]+'_hs_az315.tif'
    else:
        out_fn = os.path.splitext(fn)[0]+'_%s.tif' % product
    if not os.path.exists(out_fn):
        cmd = ['gdaldem', product]
        cmd.extend(opts)
        cmd.extend(iolib.gdal_opt_co)
        cmd.extend([fn, out_fn])
        if verbose:
            print(' '.join(cmd))
            cmd_opt = {}
        else:
            fnull = open(os.devnull, 'w')
            cmd_opt = {'stdout':fnull, 'stderr':subprocess.STDOUT}
        subprocess.call(cmd, shell=False, **cmd_opt)

    if returnma:
        ds = gdal.Open(out_fn, gdal.GA_ReadOnly)
        bma = iolib.ds_getma(ds, 1)
        return bma 
    else:
        return out_fn