Python geometry.shape函数代码示例

def main(infile, outfile, driver):

    with fio.open(infile) as src:
        meta = src.meta
        meta['driver'] = driver

    with fio.open(infile) as src, fio.open(outfile, 'w', **meta) as dst:
        with click.progressbar(src) as features:
            for feat in features:

                east = deepcopy(feat)
                west = deepcopy(feat)

                east_geom = shape(east['geometry'])
                west_geom = shape(west['geometry'])

                # if 'Point' not in asShape(feat['geometry']).type:
                #     east_geom = east_geom.simplify(0.0001).buffer(0)
                #     west_geom = west_geom.simplify(0.0001).buffer(0)

                east_geom = translate(east_geom, xoff=180)
                west_geom = translate(west_geom, xoff=-180)

                if not east_geom.is_empty:
                    east['geometry'] = mapping(east_geom)
                    dst.write(east)

                if not west_geom.is_empty:
                    west['geometry'] = mapping(west_geom)
                    dst.write(west)

    def test_PointInPolygon(self):
        """
        Test that a feature's point geometry is within its polygon geometry.
        NOTE: A feature that does not contain a polygon geometry will be skipped.
        """
        features_point_outside_polygon = {}
        for record in self.test_data:
            geoms = record['geometry']['geometries']
            pt_geom = None
            poly_geom = None
            # Check for presence of point and polygon geoms.
            for geom in geoms:
                if geom['type'] == 'Point':
                    pt_geom = shape(geom)
                if geom['type'] in ('MultiPolygon', 'Polygon'):
                    poly_geom = shape(geom)
            if pt_geom and poly_geom:
                # Check if point is within polygon. If not, flag it.
                try:
                    if not pt_geom.within(poly_geom):
                        composite_key = self._get_comp_primary_key(record)
                        features_point_outside_polygon[composite_key] = geoms
                except Exception as e:
                    continue

        json_data = json.dumps(features_point_outside_polygon)
        return json_data

def zips_in_DMA(one_dma_id, to_txt = False):

    #for both zip codes and DMAs, extract properties and geometry from each geojson feature, add each one as a dictionary entry with the DMA ID or zip # as the key and a shapely object as the value
    zips_ = {}
    for i in zips.features:
        zips_[str(i['properties']['ZCTA5CE10'])] = shape(geojson.MultiPolygon(i['geometry']['coordinates']))

    dmas_errors = []
    dmas_ = {}
    for i in nielsen.features:
        try:
            dmas_[str(i['properties']['id'])] = shape(geojson.Polygon(i['geometry']['coordinates']))
        except:
            try:
                dmas_[str(i['properties']['id'])] = shape(geojson.Polygon(i['geometry']['coordinates'][0]))
            except:
                print "error"

    dmaziplist = []
    for i in zips_.iteritems():
        if dmas_[one_dma_id].intersects(i[1]):
            dmaziplist.append(i[0])
        else:
            pass

    dma_filename = dmametadata.ix[one_dma_id]['dma_name'].replace(',', '').replace(' ', '') + one_dma_id
    if to_txt:
        with open(dma_filename + '.txt', 'wb') as zipfile:
            for i in dmaziplist:
                zipfile.write(i)
                zipfile.write('\n')

    return dmaziplist

def read_geofabric_data(netGDB,bbox=None):
    catch = {}
    with fiona.open(netGDB, layer='AHGFCatchment') as c:
        for feat in c.items(bbox=bbox):
            geom = shape(feat[1]['geometry'])
            cid = feat[1]['properties']['HydroID']
            assert cid not in catch #shouldnt be duplicates 
            catch[cid] = geom
    
    DG=nx.MultiDiGraph()
    with fiona.open(netGDB, layer='AHGFNetworkStream') as c:
        for feat in c.items(bbox=bbox):
            streamLink = shape(feat[1]['geometry'])
             #for some reason these are coming in as multipart features with only one part - no need for this
            assert streamLink.type == 'MultiLineString'
            assert len(streamLink.geoms) == 1
            streamLink = streamLink.geoms[0]
            
            ##remove this - just here for testing
            #if streamLink.representative_point().y > -40.5: #tasmania for testing
            #    continue
            
            sid = feat[1]['properties']['HydroID']
            cid = feat[1]['properties']['DrainID']
            fid = feat[1]['properties']['From_Node']
            tid = feat[1]['properties']['To_Node']
            subCatch = catch.get(cid,Polygon())
            DG.add_edge(fid, tid, id=sid,cid=cid,subCatch=subCatch,stream=streamLink)
            
    return DG

    def apply_shapely(self, method, args=None, call=True, out_geomtype=None,
                      **kwargs):
        coll = self.collection()
        out_schema = coll.schema.copy()
        if not args:
            args = []
        if out_geomtype:
            out_schema['geometry'] = out_geomtype

        tempds = self.tempds(method)
        with fiona.collection(tempds, "w", "ESRI Shapefile",
                              out_schema, crs=self.crs) as out_collection:
            for in_feature in coll:
                out_feature = in_feature.copy()
                if call:
                    geom = mapping(
                        getattr(shape(in_feature['geometry']),
                                method)(*args, **kwargs)
                    )
                else:
                    # it's not a method, it's a property
                    geom = mapping(
                        getattr(shape(in_feature['geometry']), method)
                    )

                out_feature['geometry'] = geom
                out_collection.write(out_feature)
        return Layer(tempds)

    def import_shoreline_file(self, lake_name, shoreline_file):
        """ Load the shoreline from GIS file.

        NB: Currently has side effects, loading crs and properties traits.
        """
        with fiona.open(shoreline_file) as f:
            crs = f.crs
            geometries = []
            for rec in f:
                geometries.append(rec['geometry'])
            # XXX: assuming that the properties aren't varying by geometry
            properties = rec['properties']

            if len(geometries) == 1:
                geom = shape(geometries[0])
            else:
                # XXX: this assumes we'll always get lines, not polygons or other
                geom = MultiLineString([
                    shape(geometry) for geometry in geometries])

        with self._open_file('a') as f:
            shoreline_group = self._get_shoreline_group(f)
            shoreline_group._v_attrs.crs = self._safe_serialize(crs)
            shoreline_group._v_attrs.lake_name = self._safe_serialize(lake_name)
            shoreline_group._v_attrs.original_shapefile = self._safe_serialize(shoreline_file)
            shoreline_group._v_attrs.properties = self._safe_serialize(properties)
            geometry_str = self._safe_serialize(mapping(geom))

            self._write_array(f, shoreline_group, 'geometry', np.array(geometry_str))
            f.flush()

def process_web_inputs(task, inputs):
    for param in task.inputs:
        if param.name in inputs:
            if isinstance(param, (params.RasterParameter, params.NdArrayParameter)):
                inputs[param.name] = params.RegisteredDatasetParameter(param.name).clean(inputs[param.name])

            elif isinstance(param, params.FeatureParameter):
                try:
                    inputs[param.name] = shape(inputs[param.name])
                except (ValueError, AttributeError, KeyError):
                    raise ParameterNotValidError

            elif isinstance(param, params.ListParameter):
                if isinstance(param.param_type, (params.RasterParameter, params.NdArrayParameter)):
                    inputs[param.name] = [
                        params.RegisteredDatasetParameter(param.name).clean(x) for x in inputs[param.name]
                    ]

                elif isinstance(param.param_type, params.FeatureParameter):
                    try:
                        inputs[param.name] = [shape(x) for x in inputs[param.name]]
                    except (ValueError, AttributeError, KeyError):
                        raise ParameterNotValidError

    return task.validate_inputs(inputs)

    def test_PointInPolygonGeom(self):
        """
        Test that a feature's point geometry is within its polygon geometry.
        NOTE: A feature that does not contain a polygon geometry will be skipped.
        """
        features_point_outside_polygon = []
        for record in self.test_data:
            geoms = record['geometry']['geometries']
            pt_geom = None
            poly_geom = None
            # Check for presence of point and polygon geoms.
            for geom in geoms:
                if geom['type'] == 'Point':
                    pt_geom = shape(geom)
                if geom['type'] in ('MultiPolygon', 'Polygon'):
                    poly_geom = shape(geom)
            if pt_geom and poly_geom:
                # Check if point is within polygon. If not, flag it.
                try:
                    if not pt_geom.within(poly_geom):
                        composite_key = self._get_comp_primary_key(record)
                        features_point_outside_polygon.append(composite_key)
                except Exception as e:
                    continue

        # Raise error if we don't have an empty error list
        self.assertEqual(len(features_point_outside_polygon), 0,
                         msg="The following features have point geoms outside of polygon geoms %s" % features_point_outside_polygon)

def build_aggregate(level, code, name, territories):
    print 'Building aggregate "{0}" (level={1}, code={2})'.format(name, level, code)
    polygons = []
    for tlevel, tcode in territories:
        try:
            territory = Territory.objects.get(level=tlevel, code=tcode)
        except Territory.DoesNotExist:
            print 'Territory {0}/{1} not found'.format(tlevel, tcode)
            continue
        if not shape(territory.geom).is_valid:
            print 'Skipping invalid polygon for {0}'.format(territory.name)
            continue
        if shape(territory.geom).is_empty:
            print 'Skipping empty polygon for {0}'.format(territory.name)
            continue
        polygons.append(territory.geom)

    geom = cascaded_union([shape(p) for p in polygons])
    if geom.geom_type == 'Polygon':
        geom = MultiPolygon([geom])
    try:
        territory = Territory.objects.get(level=level, code=code)
    except Territory.DoesNotExist:
        territory = Territory(level=level, code=code)
    territory.name = name
    territory.geom = geom.__geo_interface__
    territory.save()

    def test_with_initial(self):
        Fake, FakeForm = self.factory()
        territories = [TerritoryFactory() for _ in range(3)]

        fake = Fake(spatial=SpatialCoverage(
            territories=[t.reference() for t in territories],
            granularity=SPATIAL_GRANULARITIES.keys()[1]
        ))

        territory = TerritoryFactory()
        data = MultiDict({
            'spatial-territories': str(territory.id),
            'spatial-granularity': VALID_GRANULARITY
        })

        form = FakeForm(data, fake)

        form.validate()
        self.assertEqual(form.errors, {})

        form.populate_obj(fake)

        self.assertEqual(len(fake.spatial.territories), 1)
        self.assertEqual(fake.spatial.territories[0], territory.reference())
        self.assertTrue(shape(fake.spatial.geom).equals(shape(territory.geom)))

def segQuality(inVector, inImage):
    # open the vector
    lyr = fiona.open(inVector)
    features = [x for x in lyr]
    values = np.zeros([len(features), 5], dtype=float)
    # loop over features
    for  i in range(len(features)):
        geometry1 = shape(features[i]['geometry'])
        
        restFeatures = features[:i] + features[(i+ 1):]
        
        value = zonal_stats(geometry1, inImage, stats=['count'], add_stats={'mymean':mymean, "myvarianza":varianza } )
            
        df =  pd.DataFrame.from_dict(value, orient='columns', dtype=None)
        
            
        for j in range(len(restFeatures)):        
            geometry2 = shape(features[j]['geometry'])   
            if geometry2.intersects(geometry1) == True:
                #print("They touch")
                value = zonal_stats(geometry2, inImage, stats=['count'], add_stats={'mymean':mymean, "myvarianza":varianza } )
                df = df.append(pd.DataFrame.from_dict(value, orient='columns', dtype=None))
            
        values[i,0] = df.iloc[0,0] # count
        values[i,1] = df.iloc[0,1] # mean
        values[i,2] = df.iloc[0,2] # myvarianza
        values[i,3] = np.var(df.iloc[:,1]) # varianza between
        values[i,4] = len(df.iloc[1:]) # neighbours
    # get overal values   
    intraVarWeighted = np.nansum( values[:,0]*values[:,2] ) / np.nansum(values[:,0])
    interVarWeighted = np.nansum( values[:,4]*values[:,3] ) / np.nansum(values[:,4])
    normVariance = (intraVarWeighted - interVarWeighted) / (intraVarWeighted + interVarWeighted)
    numberSegments =  len(values[:,4])
    
    return( intraVarWeighted, interVarWeighted, normVariance, numberSegments )

def makeShape(code): #takes a 2-letter code
   geom = findCountry(code)
   if geom['type'] == 'Polygon':
      p=Polygon(shape(geom))
   elif geom['type'] == 'MultiPolygon':
      p=MultiPolygon(shape(geom))
   return p

def get_light_rail_sections():
    """"""

    # units derived from srs of input data
    buffer_dist = 500

    sect_dict = dict()
    sections = fiona.open(SECTIONS)
    sect_ix = generate_spatial_index(sections)

    with fiona.open(RAIL_LINES) as rail_lines:
        for rail_id, rail_feat in rail_lines.items():
            rail_type = rail_feat['properties']['TYPE']

            if 'MAX' in rail_type:
                rail_geom = shape(rail_feat['geometry'])
                rail_buffer = rail_geom.buffer(buffer_dist)

                for sect_id in sect_ix.intersection(rail_buffer.bounds):
                    # for some reason fiona sees integers and long
                    # integers of the same value as not equal
                    sect_id = int(sect_id)

                    if sect_id not in sect_dict:
                        sect_feat = sections[sect_id]
                        sect_geom = shape(sect_feat['geometry'])

                        if rail_buffer.intersects(sect_geom):
                            sect_name = sect_feat['properties']['SECTION']
                            sect_dict[sect_id] = sect_name

    return sect_dict.values()

def location_listener(self, name, location):
    global telem_data, last_point, vehicle, last_location
    telem_data['params']['location'] = {"lat": location.global_frame.lat,
                              "lon": location.global_frame.lon,
                              "alt": location.global_frame.alt}
    # print vehicle.location.global_frame
    current_point = Point(location.global_frame.lon, location.global_frame.lat)
    telem_data['params']['nofly_zone_status']['inside'] = False
    telem_data['params']['nofly_zone_status']['indexes'] = []
    for i, nofly_zone in enumerate(nofly_polygons):
        poly = shape(nofly_zone['data']['geometry'])
        if poly.contains(current_point):
            telem_data['params']['nofly_zone_status']['inside'] = True
            telem_data['params']['nofly_zone_status']['indexes'].append(nofly_zone['id'])
    further_check = False
    try:
        if global_multipoly.contains(current_point):
            further_check = True
    except:
        further_check = True
    if further_check:
        for i, nofly_zone in enumerate(global_nofly_polygons):
            poly = shape(nofly_zone['data']['geometry'])
            if poly.contains(current_point):
                telem_data['params']['nofly_zone_status']['inside'] = True
                telem_data['params']['nofly_zone_status']['indexes'].append(nofly_zone['id'])
    if telem_data['params']['nofly_zone_status']['inside'] is False:
        last_point = current_point
        last_location = location.global_frame

def checkPos(lat, lon):
    with open('5_mile_airport.json', 'r') as f:
        js = json.load(f)
    point = Point(lon, lat)
    for feature in js['features']:
        polygon = shape(feature['geometry'])
        if polygon.contains(point):
            return json.JSONEncoder().encode({ "safe" : "false"})
    with open('unitedstates.json', 'r') as f:
        js = json.load(f)
    inUS = 0
    for feature in js['features']:
        polygon = shape(feature['geometry'])
        if polygon.contains(point):
            inUS = 1
            break;
    if inUS == 1:
        with open('us_military.json', 'r') as f:
            js = json.load(f)
        for feature in js['features']:
            polygon = shape(feature['geometry'])
            if polygon.contains(point):
                return json.JSONEncoder().encode({ "safe" : "false"})
        with open('us_national_park.json', 'r') as f:
            js = json.load(f)
        for feature in js['features']:
            polygon = shape(feature['geometry'])
            if polygon.contains(point):
                return json.JSONEncoder().encode({ "safe" : "false"})
        return json.JSONEncoder().encode({ "safe" : "true"})
    else:
        return json.JSONEncoder().encode({ "safe" : "caution"})