python类shape()的实例源码

def polygonize_sk(mask, level):
    contours = measure.find_contours(mask, level)
    polys = []
    for contour in contours:
        if contour.shape[0] < 4:
            continue
        poly = Polygon(shell=contour[:, [1, 0]])
        polys.append(poly)
    polys = MultiPolygon(polys)
    return polys

def get_state_from_coordinates(coordnates):

    #coordinates = np.array([(34, -118), (40.7, -74)])
    sf = shapefile.Reader('./data/us_states_st99/st99_d00')

    #sf = shapefile.Reader("./data/states/cb_2015_us_state_20m")
    shapes = sf.shapes()
    #shapes[i].points
    fields = sf.fields
    records = sf.records()
    state_polygons = defaultdict(list)
    for i, record in enumerate(records):
        state = record[5]
        points = shapes[i].points
        poly = shape(shapes[i])
        state_polygons[state].append(poly)

    coor_state = {}
    for i in range(coordnates.shape[0]):
        lat, lon = coordnates[i]
        for state, polies in state_polygons.iteritems():
            for poly in polies:
                point = Point(lon, lat)
                if poly.contains(point):
                    coor_state[(lat, lon)] = state.lower()
    return coor_state

def subset_s2(path, file, aoi):
    # Read Data________________________________________________________________
    print("SUBSET: Read Product...")
    sentinel = ProductIO.readProduct(path+file)
    print("SUBSET: Done reading!")
    # Get Band Names and print info
    name = sentinel.getName()
    print("SUBSET: Image ID:        %s" % name)
    band_names = sentinel.getBandNames()
    print("SUBSET: Bands:       %s" % (list(band_names)))
    # Preprocessing ___________________________________________________________
    # Resampling
    parameters = HashMap()
    parameters.put('targetResolution', 10)
    print("SUBSET: resample target resolution: 10m")
    product_resample = snappy.GPF.createProduct('Resample', parameters, sentinel)
    # Geojson to wkt
    with open(aoi) as f:
            gjson = json.load(f)       
    for feature in gjson['features']:
        polygon = (feature['geometry'])
    str_poly = json.dumps(polygon)
    gjson_poly = geojson.loads(str_poly)
    poly_shp = shape(gjson_poly)
    wkt = poly_shp.wkt
    # Subset
    geom = WKTReader().read(wkt)
    op = SubsetOp()
    op.setSourceProduct(product_resample)
    op.setGeoRegion(geom)
    product_sub = op.getTargetProduct()            
    # Write Data_______________________________________________________            
    print("SUBSET: Writing subset.")
    subset = path + name + '_subset_'
    ProductIO.writeProduct(product_sub, subset, "BEAM-DIMAP")    
    print("SUBSET: Done and saved in %s" % path)

def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
    """
    Take a sequence of geographic shapes and project them to a bounds-centered
    orthographic projection.

    Args:
        shapes: Sequence of shapes, as read in by fiona.collection().
        xmin: Eastern boundary of all shapes.
        xmax: Western boundary of all shapes.
        ymin: Southern boundary of all shapes.
        ymax: Northern boundary of all shapes.

    Returns:
       Tuple of
           - Input sequence of shapes, projected to orthographic
           - PyProj projection object used to transform input shapes
    """
    latmiddle = ymin + (ymax - ymin) / 2.0
    lonmiddle = xmin + (xmax - xmin) / 2.0
    projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
               '+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
    proj = pyproj.Proj(projparams=projstr)
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)

    pshapes = []
    for tshape in shapes:
        if tshape['geometry']['type'] == 'Polygon':
            pshapegeo = shape(tshape['geometry'])
        else:
            pshapegeo = shape(tshape['geometry'])
        pshape = transform(project, pshapegeo)
        pshapes.append(pshape)  # assuming here that these are simple polygons

    return (pshapes, proj)

def sampleNoPoints(sampleimg, N, xvar, yvar):
    '''
    Sample from our "no" sample grid, where that grid contains 1s where no
    pixels should be sampled from, and 0s where they should not.

    Args:
        sampleimg: Grid of shape (len(yvar),len(xvar)) where 1's represent
            pixels from which "no" values can be sampled.
        N: Number of pixels to sample (without replacement from sampleimg.
        xvar: Numpy array of centers of columns of sampling grid.
        yvar: Numpy array of centers of rows of sampling grid.

    Returns:
        Tuple of
          - nopoints sequence of x,y tuples representing "no" samples.
          - Modified version of input sampleimg, with nopoints pixels set to 0.

    '''

    # get N points from sampleimg without replacement
    # avoid nosampleidx indices
    # return an sequence of X,Y tuples from those indices
    npixels = len(xvar) * len(yvar)
    allidx = np.arange(0, npixels)
    sampleimg = sampleimg.flatten()  # flatten out the input image
    sampleidx = np.random.choice(allidx[sampleimg == 1], size=N, replace=False)
    sampleidx.sort()
    sampleimg[sampleidx] = 0
    sampleimg.shape = (len(yvar), len(xvar))
    sampley, samplex = np.unravel_index(sampleidx, sampleimg.shape)
    xp = xvar[samplex]
    yp = yvar[sampley]
    nopoints = list(zip(xp, yp))
    return (nopoints, sampleimg, sampleidx)

def getProjectedPatch(polygon, m, edgecolor, facecolor, lw=1., zorder=10):
    polyjson = mapping(polygon)
    tlist = []
    for sequence in polyjson['coordinates']:
        lon, lat = list(zip(*sequence))
        x, y = m(lon, lat)
        tlist.append(tuple(zip(x, y)))
    polyjson['coordinates'] = tuple(tlist)
    ppolygon = shape(polyjson)
    patch = PolygonPatch(ppolygon, facecolor=facecolor, edgecolor=edgecolor,
                         zorder=zorder, linewidth=lw, fill=True, visible=True)
    return patch

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def load_map(map_path):
    sh = shapefile.Reader(map_path)
    feature = sh.shapeRecords()[0]
    first = feature.shape.__geo_interface__
    shp_geom = shape(first)
    return shp_geom

def from_db_value(self, value, expression, connection, context):
        if value is None:
            return value
        return shape(json.loads(value))

def to_python(self, value):
        if value is None or value == '':
            return None
        try:
            geometry = shape(json.loads(value))
        except Exception:
            raise ValidationError(_('Invalid GeoJSON.'))
        self._validate_geomtype(geometry)
        try:
            geometry = clean_geometry(geometry)
        except Exception:
            raise ValidationError(_('Could not clean geometry.'))
        self._validate_geomtype(geometry)
        return geometry

def calc_area_of_polygons(clons, clats):
    """
    Calculate the area of lat-lon polygons.

    We project sphere onto a flat surface using an equal area projection
    and then calculate the area of flat polygon.

    This is slow we should do some caching to avoid recomputing.
    """

    areas = np.zeros(clons.shape[1:])
    areas[:] = np.NAN

    for j in range(areas.shape[0]):
        for i in range(areas.shape[1]):

            lats = clats[:, j, i]
            lons = clons[:, j, i]

            lat_centre = lats[0] + abs(lats[2] - lats[1]) / 2
            lon_centre = lons[0] + abs(lons[1] - lons[0]) / 2

            pa = pyproj.Proj(proj_str.format(lat_centre, lon_centre))
            x, y = pa(lons, lats)

            cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
            areas[j, i] = shape(cop).area

    assert(np.sum(areas) is not np.NAN)
    assert(np.min(areas) > 0)

    return areas

def add_poly_patches(geo_df, basemap):
        geo_df['shape'] = Choroplether.make_shapes(geo_df, basemap)
        geo_df['patches'] = Choroplether.make_default_polygon_patches(
            geo_df['shape']
        )
        return geo_df

def make_shapes(geo_df, bmap):
        return geo_df['geometry'].map(lambda x: transform(bmap, shape(x)))

def rasterize_geom(geom, shape, affine, all_touched=False):
    """
    Parameters
    ----------
    geom: GeoJSON geometry
    shape: desired shape
    affine: desired transform
    all_touched: rasterization strategy

    Returns
    -------
    ndarray: boolean
    """
    geoms = [(geom, 1)]
    rv_array = features.rasterize(
        geoms,
        out_shape=shape,
        transform=affine,
        fill=0,
        dtype='uint8',
        all_touched=all_touched)

    return rv_array.astype(bool)


# https://stackoverflow.com/questions/8090229/
#   resize-with-averaging-or-rebin-a-numpy-2d-array/8090605#8090605

def rebin_sum(a, shape, dtype):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).sum(-1, dtype=dtype).sum(1, dtype=dtype)

def rasterize_pctcover_geom(geom, shape, affine, scale=None, all_touched=False):
    """
    Parameters
    ----------
    geom: GeoJSON geometry
    shape: desired shape
    affine: desired transform
    scale: scale at which to generate percent cover estimate

    Returns
    -------
    ndarray: float32
    """
    if scale is None:
        scale = 10

    min_dtype = min_scalar_type(scale ** 2)

    new_affine = Affine(affine[0]/scale, 0, affine[2],
                        0, affine[4]/scale, affine[5])

    new_shape = (shape[0] * scale, shape[1] * scale)

    rv_array = rasterize_geom(geom, new_shape, new_affine, all_touched=all_touched)

    rv_array = rebin_sum(rv_array, shape, min_dtype)

    return rv_array.astype('float32') / (scale**2)

def iter_latlong(self, indices=None):
        """Iterate over dataset as Shapely geometries in WGS 84 CRS."""
        _ = partial(project, from_proj=self.crs, to_proj='')
        if indices is None:
            for index, feature in enumerate(self):
                yield (index, _(shape(feature['geometry'])))
        else:
            for index in indices:
                yield (index, _(shape(self[index]['geometry'])))