code: def get_icemask(ds, glac_shp_fn=None):
      
      print("Masking glaciers")
      if glac_shp_fn is None:
          glac_shp_fn = get_glacier_poly()
  
      if not os.path.exists(glac_shp_fn):
          print("Unable to locate glacier shp: %s" % glac_shp_fn)
      else:
          print("Found glacier shp: %s" % glac_shp_fn)
  
      #All of the proj, extent, handling should now occur in shp2array
      icemask = geolib.shp2array(glac_shp_fn, ds)
      return icemask
docs: Generate glacier polygon raster mask for input Dataset res/extent