裁剪（求异）
#!/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 )