Reprojections




Pay Notebook Creator: Salah Ahmed0
Set Container: Numerical CPU with TINY Memory for 10 Minutes 0
Total0
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 [4]:
from osgeo import ogr, osr
from shapely import wkt
from shapely.geometry import GeometryCollection
import os
In [5]:
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 [6]:
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]:
0
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]:
9
In [28]:
GeometryCollection(lines + [wkt.loads(inters.ExportToWkt())])
Out[28]:
<shapely.geometry.collection.GeometryCollection at 0x7fe9c6c4cc18>
In [29]:
del navajo_ds
del ds