裁剪（求异）
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import json
from os.path import realpath
from shapely.geometry import MultiPolygon
from shapely.geometry import asShape
from shapely.wkt import dumps
# define our files input and output locations
input_fairways = realpath("../geodata/pebble-beach-fairways-3857.geojson")
input_greens = realpath("../geodata/pebble-beach-greens-3857.geojson")
output_wkt_sym_diff = realpath("ol3/data/ch06-01_results_sym_diff.js")
# open and load our geojson files as python dictionary
with open(input_fairways) as fairways:
fairways_data = json.load(fairways)
with open(input_greens) as greens:
greens_data = json.load(greens)
# create storage list for our new shapely objects
fairways_multiply = []
green_multply = []
# 将GeoJSON中的geometry转化为shapely中的geometry
# create shapely geometry objects for fairways
for feature in fairways_data[‘features‘]:
shape = asShape(feature[‘geometry‘])
fairways_multiply.append(shape)
# 将GeoJSON中的geometry转化为shapely中的geometry
# create shapely geometry objects for greens
for green in greens_data[‘features‘]:
green_shape = asShape(green[‘geometry‘])
green_multply.append(green_shape)
#创建MultiPolygon
# create shapely MultiPolygon objects for input analysis
fairway_plys = MultiPolygon(fairways_multiply)
greens_plys = MultiPolygon(green_multply)
#执行裁剪操作
# run the symmetric difference function creating a new Multipolygon
result = fairway_plys.symmetric_difference(greens_plys)
#将要素转化为WKT格式输出
# write the results out to well known text (wkt) with shapely dump
def write_wkt(filepath, features):
with open(filepath, "w") as f:
# create a js variable called ply_data used in html
# Shapely dumps geometry out to WKT
f.write("var ply_data = ‘" + dumps(features) + "‘")
# write to our output js file the new polygon as wkt
write_wkt(output_wkt_sym_diff, result)
合并（无dissolve）
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import json
from os.path import realpath
import shapefile # pyshp
from geojson import Feature, FeatureCollection
from shapely.geometry import asShape, MultiPolygon
from shapely.ops import polygonize
from shapely.wkt import dumps
def create_shapes(shapefile_path):
"""
Convert Shapefile Geometry to Shapely MultiPolygon
:param shapefile_path: path to a shapefile on disk
:return: shapely MultiPolygon
"""
in_ply = shapefile.Reader(shapefile_path)
# using pyshp reading geometry
ply_shp = in_ply.shapes()
ply_records = in_ply.records()
ply_fields = in_ply.fields
print(ply_records)
print(ply_fields)
if len(ply_shp) > 1:
#将OGR中的Geometry转化为shapely中的Geometry
# using python list comprehension syntax
# shapely asShape to convert to shapely geom
ply_list = [asShape(feature) for feature in ply_shp]
#转化为MultiPolygon
# create new shapely multipolygon
out_multi_ply = MultiPolygon(ply_list)
# # equivalent to the 2 lines above without using list comprehension
# new_feature_list = []
# for feature in features:
# temp = asShape(feature)
# new_feature_list.append(temp)
# out_multi_ply = MultiPolygon(new_feature_list)
print("converting to MultiPolygon: " + str(out_multi_ply))
else:
print("one or no features found")
shply_ply = asShape(ply_shp)
out_multi_ply = MultiPolygon(shply_ply)
return out_multi_ply
def create_union(in_ply1, in_ply2, result_geojson):
"""
Create union polygon
:param in_ply1: first input shapely polygon
:param in_ply2: second input shapely polygon
:param result_geojson: output geojson file including full file path
:return: shapely MultiPolygon
"""
#将外环进行合并
# union the polygon outer linestrings together
outer_bndry = in_ply1.boundary.union(in_ply2.boundary)
#将外环重新生成多边形
# rebuild linestrings into polygons
output_poly_list = polygonize(outer_bndry)
out_geojson = dict(type=‘FeatureCollection‘, features=[])
# generate geojson file output
for (index_num, ply) in enumerate(output_poly_list):
#将shapely中的Geometry转化为JSON中的Geometry
feature = dict(type=‘Feature‘, properties=dict(id=index_num))
feature[‘geometry‘] = ply.__geo_interface__
out_geojson[‘features‘].append(feature)
# create geojson file on disk
json.dump(out_geojson, open(result_geojson, ‘w‘))
# create shapely MultiPolygon
ply_list = []
for fp in polygonize(outer_bndry):
ply_list.append(fp)
out_multi_ply = MultiPolygon(ply_list)
return out_multi_ply
def write_wkt(filepath, features):
"""
:param filepath: output path for new javascript file
:param features: shapely geometry features
:return:
"""
with open(filepath, "w") as f:
# create a javascript variable called ply_data used in html
# Shapely dumps geometry out to WKT
f.write("var ply_data = ‘" + dumps(features) + "‘")
def output_geojson_fc(shply_features, outpath):
"""
Create valid GeoJSON python dictionary
:param shply_features: shapely geometries
:param outpath:
:return: GeoJSON FeatureCollection File
"""
new_geojson = []
for feature in shply_features:
feature_geom_geojson = feature.__geo_interface__
myfeat = Feature(geometry=feature_geom_geojson,
properties={‘name‘: "mojo"})
new_geojson.append(myfeat)
out_feat_collect = FeatureCollection(new_geojson)
with open(outpath, "w") as f:
f.write(json.dumps(out_feat_collect))
if __name__ == "__main__":
# define our inputs
shp1 = realpath("../geodata/temp1-ply.shp")
shp2 = realpath("../geodata/temp2-ply.shp")
# define outputs
out_geojson_file = realpath("../geodata/res_union.geojson")
output_union = realpath("../geodata/output_union.geojson")
out_wkt_js = realpath("ol3/data/results_union.js")
# create our shapely multipolygons for geoprocessing
in_ply_1_shape = create_shapes(shp1)
in_ply_2_shape = create_shapes(shp2)
# run generate union function
result_union = create_union(in_ply_1_shape, in_ply_2_shape, out_geojson_file)
# write to our output js file the new polygon as wkt
write_wkt(out_wkt_js, result_union)
# write the results out to well known text (wkt) with shapely dump
geojson_fc = output_geojson_fc(result_union, output_union)
合并（有dissolve）
utils.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from math import sqrt
import shapefile
from shapely.geometry import asShape, MultiPolygon
import json
from shapely.wkt import dumps
# calculate the size of our matplotlib output
GM = (sqrt(5) - 1.0) / 2.0
W = 8.0
H = W * GM
SIZE = (W, H)
# colors for our plots as hex
GRAY = ‘#B2B3B7‘
BLUE = ‘#6699cc‘
YELLOW = ‘#ffe680‘
RED = ‘#FF1813‘
GREEN = ‘#24CD17‘
# functions slightly modified from Sean Gilles http://toblerity.org/shapely/
# used for drawing our results using matplotlib
# matplotlib makers http://matplotlib.org/api/markers_api.html
def plot_coords_line(axis, object, color=‘#00b700‘,
symbol=‘o‘, label="text", mew=1, ms=7):
"""
mew = marker edge width in points
ms = marke size in points
"""
x, y = object.xy
axis.plot(x, y, symbol, label=label, color=color,
mew=mew, ms=ms, zorder=1)
def plot_coords_lines(axis, object, color=‘#999999‘):
for linestring in object:
x, y = linestring.xy
axis.plot(x, y, ‘o‘, color=color, zorder=2)
def plot_line(axis, object, color=‘#00b700‘, ls=‘-‘,
linewidth=2, c=‘g‘):
"""
ls is the line style options :[ ‘-‘ | ‘--‘ | ‘-.‘ | ‘:‘ | ‘steps‘ | ...]
"""
x, y = object.xy
axis.plot(x, y, color=color, linewidth=linewidth, ls=ls, c=c, zorder=1)
def plot_lines(axis, object, color=‘#00b700‘):
for line in object:
x, y = line.xy
axis.plot(x, y, color=color, alpha=0.4, linewidth=1,
solid_capstyle=‘round‘, zorder=2)
def set_plot_bounds(object, offset=1.0):
"""
Creates the limits for x and y axis plot
:param object: input shapely geometry
:param offset: amount of space around edge of features
:return: dictionary of x-range and y-range values for
"""
bounds = object.bounds
x_min = bounds[0]
y_min = bounds[1]
x_max = bounds[2]
y_max = bounds[3]
x_range = [x_min - offset, x_max + offset]
y_range = [y_min - offset, y_max + offset]
return {‘xrange‘: x_range, ‘yrange‘: y_range}
def create_shapes(shapefile_path):
"""
Convert Shapefile Geometry to Shapely MultiPolygon
:param shapefile_path: path to a shapefile on disk
:return: shapely MultiPolygon
"""
in_ply = shapefile.Reader(shapefile_path)
# using pyshp reading geometry
ply_shp = in_ply.shapes()
# ply_shp = in_ply.shapeRecords()
ply_records = in_ply.records()
ply_fields = in_ply.fields
# print ply_records
# print ply_fields
if len(ply_shp) > 1:
# using python list comprehension syntax
# shapely asShape to convert to shapely geom
ply_list = [asShape(feature) for feature in ply_shp]
# create new shapely multipolygon
out_multi_ply = MultiPolygon(ply_list)
# # equivalent to the 2 lines above without using list comprehension
# new_feature_list = []
# for feature in features:
# temp = asShape(feature)
# new_feature_list.append(temp)
# out_multi_ply = MultiPolygon(new_feature_list)
print "converting to MultiPolygon: " + str(out_multi_ply)
else:
print "one or no features found"
shply_ply = asShape(ply_shp)
out_multi_ply = MultiPolygon(shply_ply)
return out_multi_ply
def shp_2_geojson_file(shapefile_path, out_geojson):
# open shapefile
in_ply = shapefile.Reader(shapefile_path)
# get a list of geometry and records
shp_records = in_ply.shapeRecords()
# get list of fields excluding first list object
fc_fields = in_ply.fields[1:]
# using list comprehension to create list of field names
field_names = [field_name[0] for field_name in fc_fields ]
my_fc_list = []
# run through each shape geometry and attribute
for x in shp_records:
field_attributes = dict(zip(field_names, x.record))
geom_j = x.shape.__geo_interface__
my_fc_list.append(dict(type=‘Feature‘, geometry=geom_j,
properties=field_attributes))
# write GeoJSON to a file on disk
with open(out_geojson, "w") as oj:
oj.write(json.dumps({‘type‘: ‘FeatureCollection‘,
‘features‘: my_fc_list}))
def shp2_geojson_obj(shapefile_path):
# open shapefile
in_ply = shapefile.Reader(shapefile_path)
# get a list of geometry and records
shp_records = in_ply.shapeRecords()
# get list of fields excluding first list object
fc_fields = in_ply.fields[1:]
# using list comprehension to create list of field names
field_names = [field_name[0] for field_name in fc_fields ]
my_fc_list = []
# run through each shape geometry and attribute
for x in shp_records:
field_attributes = dict(zip(field_names, x.record))
geom_j = x.shape.__geo_interface__
my_fc_list.append(dict(type=‘Feature‘, geometry=geom_j,
properties=field_attributes))
geoj_json_obj = {‘type‘: ‘FeatureCollection‘,
‘features‘: my_fc_list}
return geoj_json_obj
def out_geoj(list_geom, out_geoj_file):
out_geojson = dict(type=‘FeatureCollection‘, features=[])
# generate geojson file output
for (index_num, ply) in enumerate(list_geom):
feature = dict(type=‘Feature‘, properties=dict(id=index_num))
feature[‘geometry‘] = ply.__geo_interface__
out_geojson[‘features‘].append(feature)
# create geojson file on disk
json.dump(out_geojson, open(out_geoj_file, ‘w‘))
def write_wkt(filepath, shply_geom):
"""
:param filepath: output path for new javascript file
:param shply_geom: shapely geometry features
:return:
"""
with open(filepath, "w") as f:
# create a javascript variable called ply_data used in html
# Shapely dumps geometry out to WKT
f.write("var ply_data = ‘" + dumps(shply_geom) + "‘")
union_dissolve.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from shapely.geometry import MultiPolygon
from shapely.ops import cascaded_union
from os.path import realpath
from utils import create_shapes
from utils import out_geoj
from utils import write_wkt
def check_geom(in_geom):
"""
:param in_geom: input valid Shapely geometry objects
:return: Shapely MultiPolygon cleaned
"""
plys = []
for g in in_geom:
# if geometry is NOT valid
if not g.is_valid:
print("Oh no invalid geometry")
# clean polygon with buffer 0 distance trick
new_ply = g.buffer(0)
print("now lets make it valid")
# add new geometry to list
plys.append(new_ply)
else:
# add valid geometry to list
plys.append(g)
# convert new polygons into a new MultiPolygon
out_new_valid_multi = MultiPolygon(plys)
return out_new_valid_multi
if __name__ == "__main__":
# input NOAA Shapefile
shp = realpath("../geodata/temp-all-warn-week.shp")
# output union_dissolve results as GeoJSON
out_geojson_file = realpath("../geodata/ch06-03_union_dissolve.geojson")
out_wkt_js = realpath("ol3/data/ch06-03_results_union.js")
# input Shapefile and convert to Shapely geometries
shply_geom = create_shapes(shp)
# Check the Shapely geometries if they are valid if not fix them
new_valid_geom = check_geom(shply_geom)
#进行级联合并
# run our union with dissolve
dissolve_result = cascaded_union(new_valid_geom)
# output the resulting union dissolved polygons to GeoJSON file
out_geoj(dissolve_result, out_geojson_file)
write_wkt(out_wkt_js, dissolve_result)
identity
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from shapely.geometry import asShape, MultiPolygon
from utils import shp2_geojson_obj, out_geoj, write_wkt
from os.path import realpath
def create_polys(shp_data):
"""
:param shp_data: input GeoJSON
:return: MultiPolygon Shapely geometry
"""
plys = []
for feature in shp_data[‘features‘]:
shape = asShape(feature[‘geometry‘])
plys.append(shape)
new_multi = MultiPolygon(plys)
return new_multi
def create_out(res1, res2):
"""
:param res1: input feature
:param res2: identity feature
:return: MultiPolygon identity results
"""
identity_geoms = []
for g1 in res1:
identity_geoms.append(g1)
for g2 in res2:
identity_geoms.append(g2)
out_identity = MultiPolygon(identity_geoms)
return out_identity
if __name__ == "__main__":
# out two input test Shapefiles
shp1 = realpath("../geodata/temp1-ply.shp")
shp2 = realpath("../geodata/temp2-ply.shp")
# output resulting GeoJSON file
out_geojson_file = realpath("../geodata/result_identity.geojson")
output_wkt_identity = realpath("ol3/data/ch06-04_results_identity.js")
# convert our Shapefiles to GeoJSON
# then to python dictionaries
shp1_data = shp2_geojson_obj(shp1)
shp2_data = shp2_geojson_obj(shp2)
# transform our GeoJSON data into Shapely geom objects
shp1_polys = create_polys(shp1_data)
shp2_polys = create_polys(shp2_data)
#先计算不同处，再计算相同处
# run the difference and intersection
res_difference = shp1_polys.difference(shp2_polys)
res_intersection = shp1_polys.intersection(shp2_polys)
# combine the difference and intersection polygons into results
result_identity = create_out(res_difference, res_intersection)
# export identity results to a GeoJSON
out_geoj(result_identity, out_geojson_file)
# write out new javascript variable with wkt geometry
write_wkt(output_wkt_identity, result_identity )
时间： 2024-12-15 01:13:26