code: def gdaldem_wrapper(fn, product='hs', returnma=True, verbose=True):
      
      #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
docs: 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.