modifyGeoJSON

from osgeo import ogr
import json
from geojson import loads, dumps, Feature, FeatureCollection
from shapely.geometry import shape, Point, LineString

‘‘‘
shp_driver = ogr.GetDriverByName(‘ESRI Shapefile‘)
shp_dataset = shp_driver.Open(r‘../geodata/schools.shp‘)
shp_layer = shp_dataset.GetLayer()
shp_srs = shp_layer.GetSpatialRef()
‘‘‘

filePathNE = r‘D:/Project/JavaScript/LeafletJS/WebGISDemoAngularJS/data/ne.geojson‘
filePathRegion = r‘D:/Project/JavaScript/LeafletJS/WebGISDemoAngularJS/data/region.geojson‘

def readGeoJSONFileToGeoJSON(jsonfile):
     with open(jsonfile) as jsonFile:
         jsonStr = jsonFile.read()
         featureCollection = loads(jsonStr)
         #print(dumps(featureCollection))
         return featureCollection

def readGeoJSONFileToJSONObject(jsonfile):
     with open(jsonfile) as jsonFile:
         jsonObject = json.load(jsonFile)
         return jsonObject

#
def JSONObjectToShape(jsonObject):
     geometryList = []
     for feature in jsonObject[‘features‘]:
         #将GeoJSON中的Geometry转化成shapely(Geos)中的Geometry
         # create shapely shape from geojson
         shapeObj = shape(feature[‘geometry‘])
         geometryList.append(shapeObj)
         #feature[‘geometry‘] = None
     return geometryList

jsonObjectNE = readGeoJSONFileToJSONObject(filePathNE)
geometryNEList = JSONObjectToShape(jsonObjectNE)
jsonObjectRegion = readGeoJSONFileToJSONObject(filePathRegion)
geometryRegionList = JSONObjectToShape(jsonObjectRegion)

for indexRegion, region in enumerate(geometryRegionList):
     for indexNE, ne in enumerate(geometryNEList):
         isIntersect = region.intersects(ne)
         if isIntersect:
             featureNE = jsonObjectNE[‘features‘][indexNE]
             featureRegion = jsonObjectRegion[‘features‘][indexRegion]
             featureNE[‘properties‘][‘region‘] = featureRegion[‘properties‘][‘name‘]
             #if hasattr(featureRegion[‘properties‘], ‘count‘):
             if featureRegion[‘properties‘].get(‘count‘, None) is not None:              
                 featureRegion[‘properties‘][‘count‘] = featureRegion[‘properties‘][‘count‘] + ‘,‘ + str(featureNE[‘properties‘][‘count‘])
             else:
                 featureRegion[‘properties‘][‘count‘] = str(featureNE[‘properties‘][‘count‘])

print(json.dumps(jsonObjectRegion))

modifyGeoJSON

原文地址：https://www.cnblogs.com/gispathfinder/p/9690862.html