def get_nlcd_mask(nlcd_ds, filter='not_forest', out_fn=None):
    """Generate raster mask for specified NLCD LULC filter
    """
    print("Loading NLCD LULC")
    b = nlcd_ds.GetRasterBand(1)
    l = b.ReadAsArray()
    print("Filtering NLCD LULC with: %s" % filter)
    #Original nlcd products have nan as ndv
        #12 - ice
        #31 - rock
        #11 - open water, includes rivers
        #52 - shrub, <5 m tall, >20%
        #42 - evergreeen forest
    #Should use data dictionary here for general masking
    #Using 'rock+ice+water' preserves the most pixels, although could be problematic over areas with lakes
    if filter == 'rock':
        mask = (l==31)
    elif filter == 'rock+ice':
        mask = np.logical_or((l==31),(l==12))
    elif filter == 'rock+ice+water':
        mask = np.logical_or(np.logical_or((l==31),(l==12)),(l==11))
    elif filter == 'not_forest':
        mask = ~(np.logical_or(np.logical_or((l==41),(l==42)),(l==43)))
    elif filter == 'not_forest+not_water':
        mask = ~(np.logical_or(np.logical_or(np.logical_or((l==41),(l==42)),(l==43)),(l==11)))
    else:
        print("Invalid mask type")
        mask = None
    #Write out original data
    if out_fn is not None:
        print("Writing out %s" % out_fn)
        iolib.writeGTiff(l, out_fn, nlcd_ds)
    l = None
    return mask
Generate raster mask for specified NLCD LULC filter