-
Notifications
You must be signed in to change notification settings - Fork 11
Tutorial
WIP
Karta provides a lightweight set of tools for performing analyses of geographical data. The organization of karta is around a set of container classes for vector and raster data with builtin methods for common tasks. This tutorial provides a brief introduction to some of the main parts of karta.
Should you come across any mistakes, I would very much appreciate if you could let me know, or even better, provide a pull request on Github!
The following exmaples are shown using Python 3, however karta is supported on both Python 2.7+ and Python 3.2+.
-
vector data refers to data can is treated as a set of connected or disconnected vertices. Examples might be road networks, a set of borders, geophysical survey lines, or the path taken by a bottle floating in an ocean current. In karta, these data are classified as belonging to Point, Multipoint, Line or Polygon classes. Some questions that might be asked of vector data include
- which of these points are contained in this polygon?
- how many times and where do these lines intersect each other?
- what is the average distance travelled by a particle?
-
raster data, in contrast, is data that are typically thought of in terms of pixels or a grid of values covering a surface. Examples might be an elevation map, satellite image, or an upstream area map. Operations on raster data might include slope, aspect, and curvature calculations, up and downsampling, and interpolation.
-
Coordinate reference system refers to a system of relating measurements on a coordinate system to actual positions in space, e.g. on the curved surface of the Earth. karta includes very basic support of projected and geographical coordinates, but extending this system through pyproj is something I would like to accomplish in the future.
Let's experiment with some vector data.
import karta
from karta.vector import Point, Multipoint, Line, Polygon
The Point, Multipoint, Line, and Polygon classes can all be instantiated
by providing vertices, and optionally, associated data and metadata.
pt = Point((-123.1, 49.25))
mpt = Multipoint([(-122.93, 48.62),
(-123.10, 48.54),
(-122.90, 48.49),
(-122.81, 48.56)])
line = Line([(-124.35713, 49.31437),
(-124.37857, 49.31720),
(-124.39442, 49.31833),
(-124.40311, 49.31942),
(-124.41052, 49.32203),
(-124.41681, 49.32477),
(-124.42278, 49.32588)])
poly = Polygon([(-25.41, 67.03),
(-24.83, 62.92),
(-12.76, 63.15),
(-11.44, 66.82)])
print(pt)
print(mpt)
print(line)
print(poly)
Point((-123.1, 49.25))
<class 'karta.vector.guppy.Multipoint'>
[(-122.93, 48.62), (-123.1, 48.54), (-122.9, 48.49), (-122.81, 48.56)]
Properties:{}
<class 'karta.vector.guppy.Line'>
[(-124.35713, 49.31437), (-124.37857, 49.3172)...(-124.41681, 49.32477), (-124.42278, 49.32588)]
Properties:{}
<class 'karta.vector.guppy.Polygon'>
[(-25.41, 67.03), (-24.83, 62.92)...(-11.44, 66.82), (-25.41, 67.03)]
Properties:{}
Each geometrical object now contains a vertex/vertices in a cartesian plane.
We may be interested in determining whether our point is within our polygon
print(poly.contains(pt)) # False
pt2 = Point((-25, 65))
print(poly.contains(pt2)) # True
False
True
or whether our line crosses the polygon
print(line.intersects(poly)) # False
False
There are methods for computing the nearest vertex to an external point, or the nearest point on an edge to an external point:
pt = Point((0.0, 60.0))
print(poly.nearest_point_to(pt))
print(poly.nearest_on_boundary(pt))
Point((-12.76, 63.15))
Point((-12.301580009598124, 64.42454648846582))
The vertices of multiple vertex objects can be iterated through and sliced:
subline = line[2:-2]
for pt_ in subline:
print(pt_.vertex)
(-124.39442, 49.31833)
(-124.40311, 49.31942)
(-124.41052, 49.32203)
A slice that takes part of a polygon returns a line.
print(poly[:2])
<class 'karta.vector.guppy.Line'>
[(-25.41, 67.03), (-24.83, 62.92)]
Properties:{}
Points have a distance that calculates the distance to another point. However,
if we do
pt = Point((-123.1, 49.25))
pt2 = Point((-70.66, 41.52))
print(pt.distance(pt2))
53.00666508061746
this probably isn't what we wanted. Be default, geometries in Karta use a planar
cartesian coordinate system. If our positions are meant to be geographical
coordinates, then we can provide the crs argument to each geometry at
creation, as in
pt = Point((-123.1, 49.25), crs=karta.LONLAT)
pt2 = Point((-70.66, 41.52), crs=karta.LONLAT)
pt.distance(pt2)
4109559.587727985
which now gives the great circle distance between point on the Earth, in meters.
When the coordinate system is specified, all geometrical methods obey that coordinate system. We can use this to perform queries, such which American state capitols are within 2000 km of Mexico City?
capitols = karta.read_geojson_features("us-capitols.json")[0]
near_capitols = capitols.within_radius(Point((-99.13, 19.43), crs=karta.LONLAT), 2000e3)
for cap in near_capitols:
print(cap.data["n"])
Oklahoma City, Oklahoma, United States
Montgomery, Alabama, United States
Little Rock, Arkansas, United States
Tallahassee, Florida, United States
Baton Rouge, Louisiana, United States
Jackson, Mississippi, United States
Santa Fe, New Mexico, United States
Austin, Texas, United States
TODO: describe other geometries
By using the data keyword argument, additional data can be associated with a
vector geometry.
mp = Multipoint([(1, 1), (3, 1), (4, 3), (2, 2)],
data={"species": ["T. officianale", "C. tectorum",
"M. alba", "V. cracca"]})
The data can be a list or a dictionary of lists, and are propogated through subsequent operations.
pt = mp[2]
print(pt)
print(pt.data["species"])
Point((4, 3))
M. alba
Metadata at the geometry level rather than the point level can be provided using
the properties keyword argument, which accepts a dictionary. Derived
geometries
carry the properties of their parent geometry.
poly = Polygon([(-25.41, 67.03),
(-24.83, 62.92),
(-12.76, 63.15),
(-11.44, 66.82)],
properties={"geology": "volcanic",
"alcohol": "brennivin"})
print(poly[0:3])
<class 'karta.vector.guppy.Line'>
[(-25.41, 67.03), (-24.83, 62.92), (-12.76, 63.15)]
Properties:{'geology': 'volcanic', 'alcohol': 'brennivin'}
The get_coordinate_lists method provides lists of coordinates, which make is
easy to visualize a geometry.
import matplotlib.pyplot as plt
%matplotlib qt
plt.plot(*line.get_coordinate_lists())
[<matplotlib.lines.Line2D at 0x7f5a901383c8>]
Data can be read from several common formats, including ESRI shapefiles (through
bindings to the pyshp module), GeoJSON, GPX, and comma separated value tables.
Convenience functions are kept in the karta.vector.read namespace.
Each geometry has appropriate methods to save data:
mp.to_vtk("my_vtk.vtk")
line.to_shapefile("my_shapefile")
poly.to_geojson("my_json.json")
<karta.vector.geojson.GeoJSONWriter at 0x7f5a90158400>