In [ ]:

%%bash
wget https://www.dropbox.com/s/gzxe2ehk94rclnp/data.tar.gz?dl=0
tar xf data.tar.gz?dl=0
mv data* /tmp/

In [1]:

from osgeo import ogr, osr
from shapely import wkt
from shapely.geometry import GeometryCollection
import os

In [2]:

def reproject_layer(inLayer, outSpatialRef, target_ds, target_name):
    inSpatialRef = inLayer.GetSpatialRef()
    coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
    # create the output layer
    if target_ds.GetLayer(target_name):
        target_ds.DeleteLayer(target_name)
    outLayer = target_ds.CreateLayer(
        target_name, 
        outSpatialRef,
        geom_type=inLayer.GetGeomType())
    # add fields
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)
    # get the output layer's feature definition
    outLayerDefn = outLayer.GetLayerDefn()
    # loop through the input features
    inFeature = inLayer.GetNextFeature()
    while inFeature:
        # get the input geometry
        geom = inFeature.GetGeometryRef()
        # reproject the geometry
        geom.Transform(coordTrans)
        # create a new feature
        outFeature = ogr.Feature(outLayerDefn)
        # set the geometry and attribute
        outFeature.SetGeometry(geom)
        for i in range(0, outLayerDefn.GetFieldCount()):
            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
        # add the feature to the shapefile
        outLayer.CreateFeature(outFeature)
        # dereference the features and get the next input feature
        outFeature = None
        inFeature = inLayer.GetNextFeature()

In [3]:

def transform_geometry(a_geom, b_geom):
    coordTrans = osr.CoordinateTransformation(
        a_geom.GetSpatialReference(), b_geom.GetSpatialReference())
    a_geom.Transform(coordTrans)

In [7]:

navajo_ds = ogr.Open('/tmp/data/navajo')
navajo_lyr = navajo_ds.GetLayer('cb_2017_us_aiannh_500k')

In [8]:

utah_ds = ogr.Open('/tmp/data/utah', 1)

utah_lyr = utah_ds.GetLayer('Utah')

In [9]:

navajo_feat = navajo_lyr.GetFeature(158)
navajo_feat.NAME

Out[9]:

'Navajo Nation'

In [10]:

navajo_geom = navajo_feat.geometry().Clone()

In [11]:

utah_feat = utah_lyr.GetFeature(1)
utah_feat.State

Out[11]:

'Utah'

In [12]:

utah_geom = utah_feat.geometry().Clone()

In [13]:

inSpatialRef = navajo_geom.GetSpatialReference()
outSpatialRef = utah_geom.GetSpatialReference()
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

In [14]:

navajo_geom.Transform(coordTrans)

Out[14]:

In [15]:

navajo_geom.Intersection(utah_geom).GetArea()

Out[15]:

5183358607.807744

In [16]:

e_ds = ogr.Open('/tmp/data/electric_power_lines')
e_lyr = e_ds.GetLayer('Electric_Power_Transmission_Lines')
reproject_layer(e_lyr, outSpatialRef, utah_ds, 'Electric_Power_Transmission_Lines')

In [17]:

del e_ds
del utah_ds

In [18]:

ds = ogr.Open('/tmp/data/utah')

utah = ds.GetLayer('Utah')
e_lyr = ds.GetLayer('Electric_Power_Transmission_Lines')

In [19]:

utah_feat = utah.GetFeature(1)

In [20]:

utah_geom = utah_feat.geometry().Clone()
utah_geom

Out[20]:

<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0x7fe9c6c24690> >

In [21]:

GeometryCollection([wkt.loads(navajo_geom.ExportToWkt()), wkt.loads(utah_geom.ExportToWkt())])

Out[21]:

<shapely.geometry.collection.GeometryCollection at 0x7fe9c6c4c438>

In [22]:

inters = navajo_geom.Intersection(utah_geom)

In [23]:

wkt.loads(inters.ExportToWkt())

Out[23]:

<shapely.geometry.multipolygon.MultiPolygon at 0x7fe9c6c4c1d0>

In [24]:

e_lyr.SetSpatialFilter(inters)

In [25]:

geometries = []
for f in e_lyr:
    geometries.append(f.geometry().Clone())

In [26]:

lines = [wkt.loads(g.ExportToWkt()) for g in geometries]
GeometryCollection(lines)

Out[26]:

<shapely.geometry.collection.GeometryCollection at 0x7fe9c6c4c940>

In [27]:

len(geometries)

Out[27]:

In [28]:

GeometryCollection(lines + [wkt.loads(inters.ExportToWkt())])

Out[28]:

<shapely.geometry.collection.GeometryCollection at 0x7fe9c6c4cc18>

In [29]:

del navajo_ds
del ds

Pay Notebook Creator: Salah Ahmed	0
Set Container: Numerical CPU with TINY Memory for 10 Minutes	0
Total	0

Reprojections