code: def get_glacier_poly():
      
      #rgi_fn = os.path.join(datadir, 'rgi50/regions/rgi50_merge.shp')
      #Update to rgi60, should have this returned from get_rgi.sh
      rgi_fn = os.path.join(datadir, 'rgi60/regions/rgi60_merge.shp')
      if not os.path.exists(rgi_fn):
          cmd = ['get_rgi.sh',]
          sys.exit("Missing rgi glacier data source. If already downloaded, specify correct datadir. If not, run `%s` to download" % cmd[0])
          #subprocess.call(cmd)
      return rgi_fn
docs: Calls external shell script `get_rgi.sh` to fetch:
  
      Randolph Glacier Inventory (RGI) glacier outline shapefiles 
  
      Full RGI database: rgi50.zip is 410 MB
  
      The shell script will unzip and merge regional shp into single global shp
      
      http://www.glims.org/RGI/