telluric

Interactive vector and raster geospatial data manipulation in Python

Juan Luis Cano Rodríguez - 2018-05-24 @ PyBCN

https://github.com/satellogic/telluric-talks

telluric

Who am I?

  • Aerospace Engineer with a passion for orbits 🛰
  • Chair of the Python España non profit and co-organizer of PyCon Spain 🐍
  • Software Developer at Satellogic 🌍
  • Free Software advocate and Python enthusiast 🕮
  • Hard Rock lover 🎸

1. Quick intro to geospatial data

Geospatial data, geographic data, georeferenced data, or just geodata "are defined in the ISO/TC 211 series of standards as data and information having an implicit or explicit association with a location relative to the Earth"

-- https://en.wikipedia.org/wiki/Geographic_data_and_information

A geographic information system (GIS) is a system designed to capture, store, manipulate, analyze, manage, and present spatial or geographic data.

-- https://en.wikipedia.org/wiki/Geographic_information_system

Two forms: vector and raster

  • Vector data: Discrete, mathematical entities
  • Raster data: Rectangular grids

Vector vs Raster

Image by Víctor Olaya "Sistemas de Información Geográfica", CC-BY

Vector data

  • Entities: Points, Lines and Polygons
  • Common formats: ESRI Shapefiles, GeoJSON, TopoJSON, ...

Vector data

Image by M. W. Toews https://en.wikipedia.org/wiki/File:Simple_vector_map.svg, CC-BY

Raster data

  • Common formats: GeoTIFF, JPEG2000, SRTM3, ...

Raster data

Projections

TL;DR: A mess.

  • The Earth is approximately a sphere, which is not a developable surface (i.e. it cannot be flattened onto a plane without distortion)
  • Therefore, we need projections that attempt to represent the Earth on a plane, with various tradeoffs and limitations

Projections

2. Why yet another Existing software

The landscape is complicated and somewhat difficult to navigate at times:

Landscape

Summary of key libraries:

Python C/C++
pyproj ⇒ PROJ.4
Fiona ⇒ OGR
Shapely ⇒ GEOS
rasterio ⇒ GDAL

From the Fiona documentation:

with fiona.drivers():
    with fiona.open('tests/data/coutwildrnp.shp') as source:
        meta = source.meta
        meta['schema']['geometry'] = 'Point'
        with fiona.open('test_write.shp', 'w', **meta) as sink:
            for f in source.filter(bbox=(-107.0, 37.0, -105.0, 39.0)):
               f['geometry'] = {'type': 'Point',
                    'coordinates': f['geometry']['coordinates'][0][0]}
                sink.write(f)

3. telluric library

telluric is an open source (MIT) Python library to manage vector and raster geospatial data in an interactive and easy way.

Importing for interactive use is short:

In [1]:
import telluric as tl
from telluric.constants import WGS84_CRS, WEB_MERCATOR_CRS

Basic geometry definition using GeoVector: in the simplest case, the bounds and the projection (CRS)

In [2]:
gv1 = tl.GeoVector.from_bounds(
    xmin=0, ymin=40, xmax=1, ymax=41, crs=WGS84_CRS
)
print(gv1)

GeoVector(shape=POLYGON ((0 40, 0 41, 1 41, 1 40, 0 40)), crs=CRS({'init': 'epsg:4326'}))
In [3]:
gv1

/home/juanlu/Development/telluric/telluric/plotting.py:141: UserWarning: Plotting a limited representation of the data, use the .plot() method for further customization
  "Plotting a limited representation of the data, use the .plot() method for further customization")
Out[3]:

Using Shapely geometries is also allowed:

In [4]:
from shapely.geometry import Polygon

gv2 = tl.GeoVector(
    Polygon([(0, 40), (1, 40.1), (1, 41), (-0.5, 40.5), (0, 40)]),
    WGS84_CRS
)
gv2

/home/juanlu/Development/telluric/telluric/plotting.py:141: UserWarning: Plotting a limited representation of the data, use the .plot() method for further customization
  "Plotting a limited representation of the data, use the .plot() method for further customization")
Out[4]:

Access several properties:

In [5]:
gv1.centroid

/home/juanlu/Development/telluric/telluric/plotting.py:141: UserWarning: Plotting a limited representation of the data, use the .plot() method for further customization
  "Plotting a limited representation of the data, use the .plot() method for further customization")
Out[5]:
In [6]:
gv1.area  # Real area in square meters
Out[6]:
9422706289.175217
In [7]:
gv1.within(gv2)
Out[7]:
False
In [8]:
gv1.difference(gv2)

/home/juanlu/Development/telluric/telluric/plotting.py:141: UserWarning: Plotting a limited representation of the data, use the .plot() method for further customization
  "Plotting a limited representation of the data, use the .plot() method for further customization")
Out[8]:
  • Features: GeoVector + attributes
  • FeatureCollections: a sequence of features
In [9]:
gf1 = tl.GeoFeature(
    gv1,
    {'name': 'One feature'}
)
gf2 = tl.GeoFeature(
    gv2,
    {'name': 'Another feature'}
)
print(gf1)
print(gf2)

GeoFeature(Polygon, {'name': 'One feature'})
GeoFeature(Polygon, {'name': 'Another feature'})
In [10]:
fc = tl.FeatureCollection([gf1, gf2])
fc

/home/juanlu/Development/telluric/telluric/plotting.py:141: UserWarning: Plotting a limited representation of the data, use the .plot() method for further customization
  "Plotting a limited representation of the data, use the .plot() method for further customization")
Out[10]:

Raster data:

In [11]:
# This will only save the URL in memory
rs = tl.GeoRaster2.open(
    "https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
)

# These calls will fecth some GeoTIFF metadata
# without reading the whole image
print(rs.crs)
print(rs.footprint())
print(rs.band_names)

CRS({'init': 'epsg:32618'})
GeoVector(shape=POLYGON ((101984.9999999127 2826915, 339314.9999997905 2826915, 339314.9999998778 2611485, 101985.0000002096 2611485, 101984.9999999127 2826915)), crs=CRS({'init': 'epsg:32618'}))
[0, 1, 2]
In [12]:
rs
Out[12]:
In [13]:
rs.crop(rs.footprint().buffer(-50000))
Out[13]:
In [14]:
rs[200:300, 200:240]
Out[14]:

4. Internal uses at Satellogic

  • Visualize the captures from our satellites
  • Estimate the coverage and do area computations
  • Large scale geospatial raster and vector data ingestion and manipulation

5. Future work

  • 🐛 List of issues https://github.com/satellogic/telluric/issues
  • Better visualization (good performance, unified approach)
  • Simpler constructors for objects
  • Better handling of corner cases (specially the antimeridian)
  • More advanced geometrical operations (such as geodesic buffering)
  • ...Surprises

Questions?

Final