Python feature.ShapelyFeature方法代码示例

# 需要导入模块: from cartopy import feature [as 别名]
# 或者: from cartopy.feature import ShapelyFeature [as 别名]
def paint_clip(self):
        clip = self.kwargs.pop('clip')
        clip = _to_geom_geoseries(self.df, clip, "clip")
        if clip is not None:
            if self.projection is not None:
                xmin, xmax, ymin, ymax = self.ax.get_extent(crs=ccrs.PlateCarree())
                extent = (xmin, ymin, xmax, ymax)
                clip_geom = self._get_clip(extent, clip)
                feature = ShapelyFeature([clip_geom], ccrs.PlateCarree())
                self.ax.add_feature(feature, facecolor=(1,1,1), linewidth=0, zorder=2)
            else:
                xmin, xmax = self.ax.get_xlim()
                ymin, ymax = self.ax.get_ylim()
                extent = (xmin, ymin, xmax, ymax)
                clip_geom = self._get_clip(extent, clip)
                xmin, xmax = self.ax.get_xlim()
                ymin, ymax = self.ax.get_ylim()
                polyplot(
                    gpd.GeoSeries(clip_geom), facecolor='white', linewidth=0, zorder=2,
                    extent=extent, ax=self.ax
                ) 

# 需要导入模块: from cartopy import feature [as 别名]
# 或者: from cartopy.feature import ShapelyFeature [as 别名]
def plot_lonlat_map(ax, lon, lat, data, transform, extend=None, cmap="CN_ref", bounds=np.arange(-5,76,5),
                    cbar=True, orientation="vertical",cbar_ticks=None, cbar_ticklabels=None, **kwargs):
    """
    :param ax:cartopy.mpl.geoaxes.GeoAxesSubplot, it should get from cartopy, eg:plt.axes(projection=ccrs.PlateCarree())
    :param lon: lon mesh grid for data units: degree
    :param lat: lat mesh grid for data units: degree
    :param data: radar data ,dims like lat, lon
    :param transform: The transform argument to plotting functions tells Cartopy what coordinate system your data are defined in.
    :param extend: (min_lon, max_lon, min_lat, max_lat), Latitude and longitude range, units:degrees
    :param cmap: str or Colormap, optional, A Colormap instance or registered colormap name. to see cm.py!
    :param min_max: The colorbar range(vmin, vmax). If None, suitable min/max values are automatically chosen by min max of data!
    :param cmap_bins:bins of cmaps, int
    :param cbar: bool, if True, plot with colorbar,
    :param orientation: vertical or horizontal, it is vaild when cbar is True
    :param kwargs: kwargs: other arguments for pcolormesh!
    :return:  pcolor result
    """
    assert isinstance(ax, cartopy.mpl.geoaxes.GeoAxesSubplot), "axes is not cartopy axes!"
    if extend is None:
        min_lon = np.min(lon)
        max_lon = np.max(lon)
        min_lat = np.min(lat)
        max_lat = np.max(lat)
    else:
        min_lon, max_lon, min_lat, max_lat = extend

    ax.set_aspect("equal")
    cmaps = plt.get_cmap(cmap)
    norm = BoundaryNorm(bounds, ncolors=cmaps.N, clip=True)
    pm = ax.pcolormesh(lon, lat, data, transform=transform, cmap=cmap, norm=norm, zorder=4, **kwargs)
    ax.add_feature(cfeature.OCEAN.with_scale('50m'), zorder=0)
    ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', \
                                                edgecolor='none', facecolor="white"), zorder=1)
    ax.add_feature(cfeature.LAKES.with_scale('50m'), zorder=2)
    ax.add_feature(cfeature.RIVERS.with_scale('50m'), zorder=3)

    ax.add_feature(cfeature.ShapelyFeature(CN_shp_info.geometries(), transform, \
                                           edgecolor='k', facecolor='none'), linewidth=0.5, \
                   linestyle='-', zorder=5, alpha=0.8)
    parallels = np.arange(int(min_lat), np.ceil(max_lat) + 1, 1)
    meridians = np.arange(int(min_lon), np.ceil(max_lon) + 1, 1)
    ax.set_xticks(meridians, crs=transform)
    ax.set_yticks(parallels, crs=transform)
    lon_formatter = LongitudeFormatter()
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    if cbar:
        cb = plt.colorbar(mappable=pm, ax=ax, orientation=orientation)
        if cbar_ticks is None:
            ticks = bounds
        else:
            ticks = cbar_ticks
        cb.set_ticks(ticks)

    if cbar_ticklabels is not None:
        if orientation == "vertical":
            cb.ax.set_yticklabels(cbar_ticklabels)
        else:
            cb.ax.set_xticklabels(cbar_ticklabels)
    return pm 

# 需要导入模块: from cartopy import feature [as 别名]
# 或者: from cartopy.feature import ShapelyFeature [as 别名]
def map_detecs(catalog, prefix, minmag=-5, mindep=-50, title=''):
    """Make scatter plot of detections with magnitudes (if applicable)."""
    catalog = catalog[(catalog['mag'] >= minmag)
                      & (catalog['depth'] >= mindep)].copy()

    if len(catalog) == 0:
        print('\nCatalog contains no events deeper than %s.' % mindep)
        return

    # define map bounds
    lllat, lllon, urlat, urlon, _, _, _, clon = qcu.get_map_bounds(catalog)

    plt.figure(figsize=(12, 7))
    mplmap = plt.axes(projection=ccrs.PlateCarree(central_longitude=clon))
    #mplmap.set_extent([lllon, urlon, lllat, urlat], ccrs.PlateCarree())
    mplmap.coastlines('50m', facecolor='none')

    # if catalog has magnitude data
    if not catalog['mag'].isnull().all():
        bins = [0, 5, 6, 7, 8, 15]
        binnames = ['< 5', '5-6', '6-7', '7-8', r'$\geq$8']
        binsizes = [10, 25, 50, 100, 400]
        bincolors = ['g', 'b', 'y', 'r', 'r']
        binmarks = ['o', 'o', 'o', 'o', '*']
        catalog.loc[:, 'maggroup'] = pd.cut(catalog['mag'], bins,
                                            labels=binnames)

        for i, label in enumerate(binnames):
            mgmask = catalog['maggroup'] == label
            rcat = catalog[mgmask]
            lons, lats = list(rcat['longitude']), list(rcat['latitude'])
            if len(lons) > 0:
                mplmap.scatter(lons, lats, s=binsizes[i], marker=binmarks[i],
                               c=bincolors[i], label=binnames[i], alpha=0.8,
                               zorder=10, transform=ccrs.PlateCarree())

        plt.legend(loc='lower left', title='Magnitude')

    # if catalog does not have magnitude data
    else:
        lons, lats = list(catalog['longitude']), list(catalog['latitude'])
        mplmap.scatter(lons, lats, s=15, marker='x', c='r', zorder=10)

    shape_feature = ShapelyFeature(Reader(china_border_file).geometries(),
                                   ccrs.PlateCarree(), edgecolor='black',
                                   facecolor='none')
    shape_feature1 = ShapelyFeature(Reader(fault_file).geometries(),
                                   ccrs.PlateCarree(), edgecolor='red',
                                   facecolor='none')        
    mplmap.add_feature(shape_feature)
    mplmap.add_feature(shape_feature1)

    plt.title(title, fontsize=20)
    plt.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.05)

    if mindep != -50:
        plt.savefig('%s_morethan%sdetecs.png' % (prefix, mindep), dpi=300)
    else:
        plt.savefig('%s_mapdetecs.png' % prefix, dpi=300)

    plt.close()