Geometry operations#
In the previous sections, we have learned learned how to create new vector files, how to make selections, and how to make coordinate transformations using either the entire layer or individual geometries. Now, from your previous experiences in classes that teach GIS you know that one important feature of working with vector objects (regardless of points, lines, polygons) that you can do geometric operations with them.
In this chapter, we present code for the basic operations and visualizations for them. We assume here homogenous coordinate systems in case of operations involving two geometries. If in your personal work this is not the case, you need to think conceptually and practically how you can resolve this issue (e.g., which of the two geometries should be projected?). We start by defining two geometries:
For completeness purpose, you can see the geometries below:
# Point
json.loads(point.ExportToJson())
{'type': 'Point', 'coordinates': [0.75, 0.55, 0.0]}
# Line
json.loads(line.ExportToJson())
{'type': 'LineString', 'coordinates': [[0.2, 0.2, 0.0], [0.6, 0.6, 0.0]]}
# Square
son.loads(square.ExportToJson())
{'type': 'Polygon',
'coordinates': [[[0.5, 0.25, 0.0],
[1.0, 0.25, 0.0],
[1.0, 0.75, 0.0],
[0.5, 0.75, 0.0],
[0.5, 0.25, 0.0]]]}
# Circle
json.loads(circle.ExportToJson())
Show code cell output
{'type': 'Polygon',
'coordinates': [[[0.8, 0.5, 0.0],
[0.799396002941565, 0.519027175896969, 0.0],
[0.797586443849239, 0.537977736072125, 0.0],
[0.794578609178812, 0.556775373308123, 0.0],
[0.790384610418907, 0.575344396154324, 0.0],
[0.785021335322284, 0.593610033709546, 0.0],
[0.778510379904822, 0.611498736698098, 0.0],
[0.770877961485986, 0.628938473626751, 0.0],
[0.762154813120935, 0.645859020830141, 0.0],
[0.752376059849354, 0.662192245236679, 0.0],
[0.741581077259318, 0.677872378716392, 0.0],
[0.729813332935693, 0.692836282905962, 0.0],
[0.717120211431521, 0.707023703444634, 0.0],
[0.70355282346714, 0.72037751259726, 0.0],
[0.689165800125357, 0.732843939287527, 0.0],
[0.674017072871359, 0.744372785615101, 0.0],
[0.658167640283151, 0.754917628984854, 0.0],
[0.641681322431805, 0.764436009034275, 0.0],
[0.624624503900566, 0.772889598606356, 0.0],
[0.607065866477562, 0.780244358079532, 0.0],
[0.589076112598482, 0.786470672433222, 0.0],
[0.570727680652828, 0.791543470497063, 0.0],
[0.552094453300079, 0.795442325903662, 0.0],
[0.533251459970303, 0.798151539338376, 0.0],
[0.514274574747123, 0.799660201754902, 0.0],
[0.495240210849558, 0.799962238302163, 0.0],
[0.476225012942963, 0.799056432785583, 0.0],
[0.457305548518014, 0.79694643256428, 0.0],
[0.438557999580443, 0.793640733864434, 0.0],
[0.420057855892989, 0.789152647567983, 0.0],
[0.401879611004773, 0.7835002456144, 0.0],
[0.384096462292061, 0.776706288231374, 0.0],
[0.366780016218268, 0.768798132287401, 0.0],
[0.35, 0.759807621135332, 0.0],
[0.333823980840167, 0.749770956390431, 0.0],
[0.3183170938587, 0.73872855215925, 0.0],
[0.303541779816415, 0.726724872306277, 0.0],
[0.289557533688104, 0.713808251413659, 0.0],
[0.276420665097274, 0.700030700154887, 0.0],
[0.264184071577164, 0.685447695866181, 0.0],
[0.25289702557105, 0.670117959158831, 0.0],
[0.242604976029507, 0.654103217472022, 0.0],
[0.233349365403523, 0.637467956518223, 0.0],
[0.225167462770379, 0.620279160621984, 0.0],
[0.218092213764227, 0.602606042997701, 0.0],
[0.212152107915651, 0.584519767052429, 0.0],
[0.207371063934378, 0.566093159835962, 0.0],
[0.203768333397082, 0.547400418792005, 0.0],
[0.201358423228075, 0.528516812991255, 0.0],
[0.200151037285044, 0.50951838004942, 0.0],
[0.200151037285044, 0.49048161995058, 0.0],
[0.201358423228075, 0.471483187008745, 0.0],
[0.203768333397082, 0.452599581207995, 0.0],
[0.207371063934378, 0.433906840164038, 0.0],
[0.212152107915651, 0.415480232947571, 0.0],
[0.218092213764227, 0.397393957002299, 0.0],
[0.225167462770379, 0.379720839378016, 0.0],
[0.233349365403523, 0.362532043481777, 0.0],
[0.242604976029507, 0.345896782527978, 0.0],
[0.25289702557105, 0.329882040841169, 0.0],
[0.264184071577164, 0.314552304133818, 0.0],
[0.276420665097274, 0.299969299845113, 0.0],
[0.289557533688104, 0.286191748586341, 0.0],
[0.303541779816414, 0.273275127693723, 0.0],
[0.3183170938587, 0.26127144784075, 0.0],
[0.333823980840167, 0.250229043609569, 0.0],
[0.35, 0.240192378864668, 0.0],
[0.366780016218268, 0.231201867712599, 0.0],
[0.384096462292061, 0.223293711768626, 0.0],
[0.401879611004773, 0.216499754385599, 0.0],
[0.42005785589299, 0.210847352432017, 0.0],
[0.438557999580443, 0.206359266135566, 0.0],
[0.457305548518014, 0.20305356743572, 0.0],
[0.476225012942963, 0.200943567214417, 0.0],
[0.495240210849558, 0.200037761697837, 0.0],
[0.514274574747123, 0.200339798245098, 0.0],
[0.533251459970303, 0.201848460661624, 0.0],
[0.552094453300079, 0.204557674096338, 0.0],
[0.570727680652828, 0.208456529502938, 0.0],
[0.589076112598482, 0.213529327566778, 0.0],
[0.607065866477561, 0.219755641920468, 0.0],
[0.624624503900566, 0.227110401393645, 0.0],
[0.641681322431805, 0.235563990965725, 0.0],
[0.658167640283151, 0.245082371015146, 0.0],
[0.674017072871359, 0.255627214384899, 0.0],
[0.689165800125357, 0.267156060712473, 0.0],
[0.70355282346714, 0.27962248740274, 0.0],
[0.717120211431521, 0.292976296555366, 0.0],
[0.729813332935693, 0.307163717094038, 0.0],
[0.741581077259318, 0.322127621283608, 0.0],
[0.752376059849354, 0.337807754763321, 0.0],
[0.762154813120935, 0.354140979169859, 0.0],
[0.770877961485986, 0.371061526373248, 0.0],
[0.778510379904822, 0.388501263301902, 0.0],
[0.785021335322284, 0.406389966290454, 0.0],
[0.790384610418907, 0.424655603845676, 0.0],
[0.794578609178812, 0.443224626691877, 0.0],
[0.797586443849239, 0.462022263927875, 0.0],
[0.799396002941565, 0.480972824103031, 0.0],
[0.8, 0.5, 0.0],
[0.8, 0.5, 0.0]]]}
We see different geometries: (1) a sqare, (2) a circle, (3) a line, and (4) a point. Some of tehm overlap, others do not. Now, we can apply operations/comparisons you have learned in your introductory GIS class. If you are uncertain of the different operations, you can can check some of the references provided at the bottom of the page. Let’s go one-by-one over the most important operations we will use in our daily work. Again: remember that we make the assumption that the geometries are in the projection. If they weren’t any of the operations below would result in an empty result. Likewise, remember that any operations are only executable at the level of geometries. When you e.g., intersect two layers in QGIS, then what the software does under the hood is to iterate over every binary combination of two geometries between the two layers. In python (or any other programing language) we will have to work with loops.
We continue now with two types of comparisons of geometries: (1) we can compare two geometries and test their relationship, (2) we can make operations using two geometries and create a new geometry
Relationship between geometries#
Within ogr we can ask the following queries to assess the overlap between two geometries. The return value here is always True or False depending on the condition. There is a comprehensive collection of operations in the book by Garrard (page 131), here I show the example for the functions Intersect(), Within() and Contains(). Compare the outcomes to the figure above.
circle.Intersect(square)
True
square.Intersect(point)
False
line.Within(circle)
False
point.Within(circle)
True
circle.Contains(line)
False
circle.Contains(point)
True
One can use these operations to make simple tests. For example, sometimes it is only interesting to ask whether an electricity line goes through a settlement (polygon), but we don’t want to create actually a new geometry. In these cases we would would want to use the tests for relationships between geometries. Can you think of other examples?
Overlay operations#
Using overlay operations we do not ask anymore solely for the relationship between two geometries. instead, we are actually creating new geometries. In your introductory GIS class and in other courses you probably have already gained some experience with these types of operations in software environments, such as QGIS or ArcGIS. Here, we take some of these operations and
Difference#
This operation creates the difference of two geometries. The syntax here (and in other overlay operations) is always geom1.Difference(geom2). If you look at the two polygons (i…e, the circle and the square), you can see that you can build two ways of differences. One, in which you erase the circle from the square, and one in which you erase the square from the circle. Belo you find an example for both. Can you think of an example where this might be helpful/problematic?
square_diff_circle = square.Difference(circle)
circle_diff_square = circle.Difference(square)
Intersection#
This operation is probably the one you know best. Here, we basically are interested in the overlap of two geometries - and it is not important which geometry is the first and the second in the syntax of the operation:
square_int_circle = square.Intersection(circle)
Union#
The union of two geometries is what you got to know in other contexts as dissolve. Again, the order of the arguments is not important.
square_union_circle = square.Union(circle)
Other operations#
Last but not least, there are a lot of other operations, that we often use. Again, have a look into the book by Garrard for a comprehensive list. Here, we present some of the most frequently used operations
Buffer#
Buffering a geometry basically means that we create a region around a geometry at a specified distance, resulting in a new geometry. The buffer can be applied expanding, which means towards the outside of the geometry (points, lines, polygons), or shrinking, which means towawrds the inside of the geometry (can only be applied to polygons). The only argument needed is the buffer distance. Here, you need to be careful, as these has to be provided in the units of the underlying coordinate system (e.g., in a UTM coordinate system we have m, in a pseudo-mercator projection we have degrees).
Calculate area and/or length#
Often, when we have done an overlay operation and have generated a new geometry, we want to assess some properties of this geometry. For example, you have a .gpkg file of protected areas and you want to calculate the average area or the sum of the areas. Likewise, you may have a .gpkg file of roads and you want to know how many kilometers of roads are in a certain region. Below are now two examples. The important element to know here is that the output value carries the unit of the underlying reference system. Thus, when you for example calculate the area of a polygon, make sure that you have the correct projection, so that you calculate the area in m^2 and not in degrees^2. For our example geometries it does not matter at all, bbecause we have not assigned any coordinate reference system to them.
circle.Area()
0.2825535620699948
line.Length()
0.565685424949238
This is all. With this information here and the ones from the previous chapters, you are prepared for the vast majority of vector operations we commonly do in Geography. If you feel something is missing - feel free to drop us a line!