Python shapely.geometry.asShape() Examples

def parse_flickr_geojson(flickr_file):
    json_data = codecs.open(flickr_file, "r", "utf-8").read()
    
    data = demjson.decode(json_data)

    features =  data['features']
    for feature in features:
        woe_id = str(feature['properties']['woe_id'])
      
        name = feature['properties']['label']
        
        if not name:
            continue
       
        feature_type = feature['properties']['place_type']
        feature_code = str(feature['properties']['place_type_id'])
        json_geometry  = feature['geometry']
        updated = "2011-01-08 00:00:00+00"  #i.e. as from http://code.flickr.com/blog/2011/01/08/flickr-shapefiles-public-dataset-2-0/
        geometry = asShape(json_geometry).wkt

        out_line = ['F', woe_id, name, feature_type, feature_code, updated, geometry ]
        
        print "\t".join(out_line)
            
    return flickr_file 

def check_file(filename):
    with collection(filename, 'r') as features:
        for feature in features:
            try:
                shape = asShape(feature['geometry'])
                if not shape.is_valid:
                    geometry = json.dumps(feature, indent=2)
                    print "Invalid geometry:\n"
                    print geometry
                    print '\n'
            except:
                print "Error parsing:\n"
                print json.dumps(feature, indent=2) 

def chunk(featureFileName, sectionFileName, pattern, key = None):

    # Load and index
    with collection(featureFileName, "r") as featureFile:
        featureIdx = index.Index()
        features = []
        for feature in featureFile:
            try:
                shape = asShape(feature['geometry'])
                features.append(feature)
                featureIdx.add(len(features) - 1, shape.bounds)
            except ValueError:
                print "Error parsing feature"
                pprint(feature)

        # Break up by sections and export
        with collection(sectionFileName, "r") as sectionFile:
            i = 0
            for section in sectionFile:
                fileName = pattern % i
                if key:
                    fileName = pattern % section['properties'][key]
                    properties = {}
                    try:
                        with collection(fileName, 'w', 'ESRI Shapefile',
                                schema = featureFile.schema,
                                crs = featureFile.crs) as output:
                            sectionShape = asShape(section['geometry'])
                            for j in featureIdx.intersection(sectionShape.bounds):
                                if asShape(features[j]['geometry']).intersects(sectionShape):
                                    properties = features[j]['properties']
                                    output.write(features[j])
                            print "Exported %s" % fileName
                            i = i + 1
                    except ValueError:
                        print "Error exporting " + fileName
                        pprint(properties)
                        pprint(featureFile.schema) 

def test_create_extent(self):

        package = factories.Dataset()

        geojson = json.loads(self.geojson_examples['point'])

        shape = asShape(geojson)
        package_extent = PackageExtent(package_id=package['id'],
                                       the_geom=WKTElement(shape.wkt,
                                                           self.db_srid))
        package_extent.save()

        assert_equals(package_extent.package_id, package['id'])
        if legacy_geoalchemy:
            assert_equals(Session.scalar(package_extent.the_geom.x),
                          geojson['coordinates'][0])
            assert_equals(Session.scalar(package_extent.the_geom.y),
                          geojson['coordinates'][1])
            assert_equals(Session.scalar(package_extent.the_geom.srid),
                          self.db_srid)
        else:
            from sqlalchemy import func
            assert_equals(
                Session.query(func.ST_X(package_extent.the_geom)).first()[0],
                geojson['coordinates'][0])
            assert_equals(
                Session.query(func.ST_Y(package_extent.the_geom)).first()[0],
                geojson['coordinates'][1])
            assert_equals(package_extent.the_geom.srid, self.db_srid) 

def _get_extent_object(self, geometry):
        if isinstance(geometry, basestring):
            geometry = json.loads(geometry)
        shape = asShape(geometry)
        return PackageExtent(package_id='xxx',
                             the_geom=WKTElement(shape.wkt, 4326)) 

def asShape(self):
        return asShape(self) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue
        borough = ""
        boros = {"1":"Manhattan", "2": "Bronx", "3":"Brooklyn", "4": "Queens", "5": "Staten Island"}
        if properties.get("BORO"):
            borough = boros[properties["BORO"]] + ", "
        
        block = ""
        if properties.get("BLOCK"):
            block = "Block " + str(properties["BLOCK"]) +", "
        name = borough + block + "Lot " + str(properties["LOT"])

        bbl = properties.get("BBL")  #stored with uris
        
        #feature code mapping
        feature_code = "ADM7"
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name  + bbl
        uri_bbl = "bbl:"+bbl
         
        timeframe = {}
       
        updated = datetime.datetime.utcnow().replace(second=0, microsecond=0).isoformat()

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "uris":[uri, uri_bbl],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }
        #print place
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue

       
        if properties["NAME"]:
            name = properties["NAME"]
        else:
            continue
        
        #feature code mapping
        feature_code = "ADM3"
        if properties["LSAD"] == "Resvn":
            feature_code = "RESV"
   
        area = properties["CENSUSAREA"]
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name + "." + properties["GEO_ID"] + "."+ feature["id"]
         
        timeframe = {}
        timeframe = {"start": "2000-01-01","start_range":0, "end": "2010-01-01", "end_range":0}
        
        updated = "2012-01-31"

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "area": area,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }
        #print place
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        #calculate centroid
        geom_obj = asShape(geometry)
        if simplify_tolerance:
            geom_obj = geom_obj.simplify(simplify_tolerance)
        
        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue
        geometry = mapping(geom_obj)

            
        if properties["FULL_NAME"]:
            name = properties["FULL_NAME"]
                    
        #feature code mapping
        feature_code = "ADM1H"
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name + "." + properties["ID"] + "."+ str(properties["VERSION"])
        
        #1766/07/02  to 1766-01-01
        timeframe = {"start": properties["START_DATE"].replace('/','-'), "start_range":0, 
                     "end": properties["END_DATE"].replace('/','-'), "end_range":0}
        
        #TODO admin? for counties?
        
        updated = "2011-10-01"
        
        
        area = properties["AREA_SQMI"]
        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "updated": updated,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[],
            "area": area
            
        }
        
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue

       
        if properties["NAME"]:
            name = properties["NAME"]
        else:
            continue
        
        #feature code mapping
        feature_code = "ADM3"
        if properties["LSAD_TRANS"] == "Reservation":
            feature_code = "RESV"
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name + "." + properties["COUSUBFP"] + "."+ feature["id"]
         
        timeframe = {}
        timeframe = {"start": "1990-01-01", "start_range":0, "end": "2000-01-01", "end_range":0}
        
        updated = "2012-01-31"

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }
        #print place
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue

        boros = {"1":"Manhattan", "2": "Bronx", "3":"Brooklyn", "4": "Queens", "5": "Staten Island"}
        name = "Block "+str(properties["BLOCK"])
        if properties.get("BORO"):
            borough = boros[properties["BORO"]]
            name = borough + ", Block "+ str(properties["BLOCK"])
                
        
        #feature code mapping
        feature_code = "ADM6"
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        
        bb = str(properties.get("BORO"))+str(properties["BLOCK"])
        
        uri = uri_name + bb +"/"+feature["id"]
        uri_bb = "bb:"+bb
        timeframe = {}
       
        updated = datetime.datetime.utcnow().replace(second=0, microsecond=0).isoformat()

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "uris":[uri,uri_bb],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }
        #print place
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)
        if simplify_tolerance:
            geom_obj = geom_obj.simplify(simplify_tolerance)
        
        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue
        geometry = mapping(geom_obj)

            
        if properties["FULL_NAME"]:
            name = properties["FULL_NAME"]
                    
        #feature code mapping
        feature_code = "ADM2H" #default code (building)
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name + "." + properties["ID"] + "."+ str(properties["VERSION"])
    
        #1766/07/02  to 1766-01-01
        timeframe = {"start": properties["START_DATE"].replace('/','-'), "start_range":0, 
                     "end": properties["END_DATE"].replace('/','-'), "end_range":0}
        
        #TODO admin? for counties?
        
        updated = "2010-02-12"
        
        
        area = properties["AREA_SQMI"]
        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "updated": updated,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[],
            "area": area
            
        }
        
        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue

        name = properties.get("NAME")

        #feature code mapping
        feature_code = "HSTS"
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        
        
        uri = uri_name + properties["LP_Number"]

        timeframe = {}
       
        updated = datetime.datetime.utcnow().replace(second=0, microsecond=0).isoformat()

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }

        dump.write(uri, place) 

def extract_shapefile(shapefile, uri_name, simplify_tolerance=None):
    
    for feature in collection(shapefile, "r"):
        
        geometry = feature["geometry"]
        properties = feature["properties"]
        
        #calculate centroid
        geom_obj = asShape(geometry)

        try:
            centroid = [geom_obj.centroid.x , geom_obj.centroid.y]    
        except AttributeError:
            print "Error: ", feature
            continue

       
        if properties["NAME"]:
            name = properties["NAME"]
        else:
            continue
        
        #feature code mapping
        feature_code = "ADM3"
       
        area = properties["AREATOT"]
                
        source = properties  #keep all fields anyhow
        
        # unique URI which internally gets converted to the place id
        # Must be unique!
        uri = uri_name + "." + properties["GEOID"] + "."+ feature["id"]
         
        timeframe = {}
        timeframe = {"start": "1980-01-01", "start_range": 0, "end": "1990-01-01", "end_range":0}
        
        updated = "2012-01-31"

        place = {
            "name":name,
            "centroid":centroid,
            "feature_code": feature_code,
            "geometry":geometry,
            "is_primary": True,
            "source": source,
            "alternate": [],
            "updated": updated,
            "area": area,
            "uris":[uri],
            "relationships": [],
            "timeframe":timeframe,
            "admin":[]

        }
        #print place
        dump.write(uri, place) 

def test_update_extent(self):

        package = factories.Dataset()

        geojson = json.loads(self.geojson_examples['point'])

        shape = asShape(geojson)
        package_extent = PackageExtent(package_id=package['id'],
                                       the_geom=WKTElement(shape.wkt,
                                                           self.db_srid))
        package_extent.save()
        if legacy_geoalchemy:
            assert_equals(
                Session.scalar(package_extent.the_geom.geometry_type),
                'ST_Point')
        else:
            from sqlalchemy import func
            assert_equals(
                Session.query(
                    func.ST_GeometryType(package_extent.the_geom)).first()[0],
                'ST_Point')

        # Update the geometry (Point -> Polygon)
        geojson = json.loads(self.geojson_examples['polygon'])

        shape = asShape(geojson)
        package_extent.the_geom = WKTElement(shape.wkt, self.db_srid)
        package_extent.save()

        assert_equals(package_extent.package_id, package['id'])
        if legacy_geoalchemy:
            assert_equals(
                Session.scalar(package_extent.the_geom.geometry_type),
                'ST_Polygon')
            assert_equals(
                Session.scalar(package_extent.the_geom.srid),
                self.db_srid)
        else:
            assert_equals(
                Session.query(
                    func.ST_GeometryType(package_extent.the_geom)).first()[0],
                'ST_Polygon')
            assert_equals(package_extent.the_geom.srid, self.db_srid) 

def save_package_extent(package_id, geometry = None, srid = None):
    '''Adds, updates or deletes the package extent geometry.

       package_id: Package unique identifier
       geometry: a Python object implementing the Python Geo Interface
                (i.e a loaded GeoJSON object)
       srid: The spatial reference in which the geometry is provided.
             If None, it defaults to the DB srid.

       Will throw ValueError if the geometry object does not provide a geo interface.

       The responsibility for calling model.Session.commit() is left to the
       caller.
    '''
    db_srid = int(config.get('ckan.spatial.srid', '4326'))


    existing_package_extent = Session.query(PackageExtent).filter(PackageExtent.package_id==package_id).first()

    if geometry:
        shape = asShape(geometry)

        if not srid:
            srid = db_srid

        package_extent = PackageExtent(package_id=package_id,
                                       the_geom=WKTElement(shape.wkt, srid))

    # Check if extent exists
    if existing_package_extent:

        # If extent exists but we received no geometry, we'll delete the existing one
        if not geometry:
            existing_package_extent.delete()
            log.debug('Deleted extent for package %s' % package_id)
        else:
            # Check if extent changed
            if not compare_geometry_fields(package_extent.the_geom, existing_package_extent.the_geom):
                # Update extent
                existing_package_extent.the_geom = package_extent.the_geom
                existing_package_extent.save()
                log.debug('Updated extent for package %s' % package_id)
            else:
                log.debug('Extent for package %s unchanged' % package_id)
    elif geometry:
        # Insert extent
        Session.add(package_extent)
        log.debug('Created new extent for package %s' % package_id)