Python——图像的白化和裁剪

白化——maskout

# -*- encoding: utf-8 -*-
import shapefile
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch


def find_region(shp_path, area_str):
    file = shapefile.Reader(shp_path)  
    opt = 0
    for k in range(len(file.records()[0])):
        i = 0
        for shape_rec in file.shapeRecords():
            if isinstance(shape_rec.record[k], bytes):
                if str(shape_rec.record[k].decode("gbk")) == area_str:
                    return k, i, shape_rec.record[k]
            else:
                if str(shape_rec.record[k]) == area_str:
                    return k, i, shape_rec.record[k]
            i = i + 1
    if opt == 0:
        print("无此区域")


def maskout_areas(originfig, ax, shp_path, area_str, proj=None):
    if area_str == ["THE_ALL"]:
        import maskout_all
        maskout_all.maskout_areas(originfig=originfig, ax=ax, shp_path=shp_path, proj=proj)
        return 0
    region = []
    for i in range(len(area_str)):
        region_k, region_i, region_str = find_region(shp_path, area_str[i])  # type: ignore
        region.append(region_str)
    sf = shapefile.Reader(shp_path)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[int(region_k)] in region:
            pts = shape_rec.shape.points
            prt = list(shape_rec.shape.parts) + [len(pts)]
            for i in range(len(prt) - 1):
                for j in range(prt[i], prt[i + 1]):
                    if proj:
                        vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
                    else:
                        vertices.append((pts[j][0], pts[j][1]))
                codes += [Path.MOVETO]
                codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
                codes += [Path.CLOSEPOLY]
            clip = Path(vertices, codes)
            clip = PathPatch(clip, transform=ax.transData)
    for contour in originfig.collections:
        contour.set_clip_path(clip)


def maskout_points(lon, lat, shp_path, area_str):
    if area_str == ["THE_ALL"]:
        import maskout_all
        point_locs = maskout_all.maskout_points(lon=lon, lat=lat, shp_path=shp_path)
        return point_locs
    point_locs = []
    for i in range(len(area_str)):
        region_k, region_i, region_str = find_region(shp_path, area_str[i])  # type: ignore
        p = Path(shapefile.Reader(shp_path).shapes()[region_i].points)
        for j in range(len(lon)):
            if p.contains_points([[lon[j], lat[j]]]):  # type: ignore
                point_locs = point_locs + [j]
    return point_locs

裁剪——shp2clip

def shp2clip(originfig, ax, shpfile, fieldVals):
    """
    This method enables you to maskout the unneccessary data
                                            outside the interest region
    :param ax:  the Axes instance
    :param shpfile:  the shape file used for clip
    :param fieldVals:  thi features attributes value list in shape file,
                    outside the region the data is to be masked out
    :return:
    """
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[0] in fieldVals:  
            pts = shape_rec.shape.points
            prt = list(shape_rec.shape.parts) + [len(pts)]
            for i in range(len(prt) - 1):
                for j in range(prt[i], prt[i + 1]):
                    vertices.append((pts[j][0], pts[j][1])) 
                codes += [Path.MOVETO]
                codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
                codes += [Path.CLOSEPOLY]
            clip = Path(vertices, codes)
            clip = PathPatch(clip, transform=ax.transData)

    for contour in originfig.collections:
        contour.set_clip_path(clip)
    return clip