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
python类shape()的实例源码
utils.py 文件源码
项目:kaggle-dstl-satellite-imagery-feature-detection
作者: u1234x1234
项目源码
文件源码
阅读 17
收藏 0
点赞 0
评论 0
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'])))