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
Image by Víctor Olaya "Sistemas de Información Geográfica", CC-BY
Image by M. W. Toews https://en.wikipedia.org/wiki/File:Simple_vector_map.svg, CC-BY
A method to locate points on earth
TL;DR: A mess.
The landscape is complicated and somewhat difficult to navigate at times:
Summary of key libraries:
Python | C/C++ | |
---|---|---|
pyproj | ⇒ | PROJ.4 |
Fiona | ⇒ | OGR |
Shapely | ⇒ | GEOS |
rasterio | ⇒ | GDAL |
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:
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)
gush_dan = tl.GeoVector.from_bounds(
xmin=34.6, ymin=32, xmax=35.0, ymax=32.3, crs=WGS84_CRS
)
print(gush_dan)
GeoVector(shape=POLYGON ((34.6 32, 34.6 32.3, 35 32.3, 35 32, 34.6 32)), crs=CRS({'init': 'epsg:4326'}))
gush_dan
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
Using Shapely geometries is also allowed:
from shapely.geometry import Polygon
south_tlv = tl.GeoVector(
Polygon([(34.75, 32.04), (34.75, 32.06), (34.77,32.06), (34.77, 32.04), ( 34.75, 32.04)]),
WGS84_CRS
)
south_tlv
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
gush_dan.centroid
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
gush_dan.area # Real area in square meters
1259014132.9260154
gush_dan.within(south_tlv)
False
south_tlv.within(gush_dan)
True
print(south_tlv.difference(gush_dan))
GeoVector(shape=GEOMETRYCOLLECTION EMPTY, crs=CRS({'init': 'epsg:4326'}))
gush_dan.difference(south_tlv)
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
GeoVector
+ attributesFeatureCollection
s: a sequence of featuresgf1 = tl.GeoFeature(
gush_dan,
{'district': 'גוש דן'}
)
gf2 = tl.GeoFeature(
south_tlv,
{'city': 'תל אביב'}
)
print(gf1)
print(gf2)
GeoFeature(Polygon, {'district': 'גוש דן'}) GeoFeature(Polygon, {'city': 'תל אביב'})
fc = tl.FeatureCollection([gf1, gf2])
fc
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
# 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.band_names)
CRS({'init': 'epsg:32618'}) [0, 1, 2]
rs.footprint() # is the bouding box
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
rs
rs.image[:,200:300, 200:240]
masked_array( data=[[[15, 20, 23, ..., 41, 38, 29], [17, 15, 23, ..., 41, 46, 46], [18, 17, 17, ..., 39, 39, 44], ..., [106, 49, 62, ..., 46, 58, 50], [255, 255, 255, ..., 38, 37, 43], [255, 229, 255, ..., 40, 31, 35]], [[94, 98, 100, ..., 129, 98, 62], [98, 96, 100, ..., 127, 127, 133], [98, 96, 94, ..., 126, 126, 129], ..., [145, 79, 116, ..., 47, 58, 51], [255, 255, 255, ..., 43, 41, 45], [255, 249, 255, ..., 41, 33, 37]], [[131, 131, 131, ..., 161, 125, 86], [133, 133, 135, ..., 164, 160, 161], [133, 132, 129, ..., 159, 160, 157], ..., [141, 77, 108, ..., 29, 34, 31], [255, 255, 255, ..., 31, 30, 31], [255, 255, 255, ..., 23, 19, 27]]], mask=[[[False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], ..., [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False]], [[False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], ..., [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False]], [[False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], ..., [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False], [False, False, False, ..., False, False, False]]], fill_value=999999, dtype=uint8)
# Crop fetch only the data required to build the raster
rs.crop(rs.footprint().buffer(-50000))
rs[200:300, 200:240] # the image is streched due to the presentation layer
# Rasterizing Vectors
# Geting OpenStreetMap data of Tel Aviv
from shapely.geometry import shape
import json
with open("telaviv.geojson") as data:
features = json.load(data)["features"]
fc = tl.FeatureCollection([tl.GeoFeature(tl.GeoVector(shape(feature["geometry"])),feature["properties"]) for feature in features])
fc
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
streets_raster = fc.rasterize(dest_resolution=1)
streets_raster
print(gush_dan)
GeoVector(shape=POLYGON ((34.6 32, 34.6 32.3, 35 32.3, 35 32, 34.6 32)), crs=CRS({'init': 'epsg:4326'}))
projected = gush_dan.reproject(tl.constants.WEB_MERCATOR_CRS)
print(projected)
GeoVector(shape=POLYGON ((3851654.381447267 3763310.627144652, 3851654.381447267 3802755.031859114, 3896182.177764576 3802755.031859114, 3896182.177764576 3763310.627144652, 3851654.381447267 3763310.627144652)), crs=CRS({'init': 'epsg:3857'}))
streets_raster.resize(dest_height=streets_raster.width*20, dest_width=streets_raster.width*10)
A Cloud Optimized GeoTIFF (COG) is a regular GeoTIFF file, aimed at being hosted on a HTTP file server, with an internal organization that enables more efficient workflows on the cloud. It does this by leveraging the ability of clients issuing HTTP GET range requests to ask for just the parts of a file they need.1
We use COGs to enable distributed computation:
from telluric.vectors import generate_tile_coordinates
list_of_regions = list(generate_tile_coordinates(rs.footprint(), num_tiles=(10,10)))
# This line is to visualize and not required
tl.FeatureCollection(tl.GeoFeature(tile, {}) for tile in list_of_regions)
/home/guyd/pEnvs/telluric-dev/local/lib/python3.5/site-packages/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")
chunk0 = rs.crop(list_of_regions[17])
chunk0
def worker(raster_filename, roi):
raster = tl.GeoRaster2.open(raster_filename) #Lazy loading of the raster
chunk0 = raster.crop(roi)
return chunk0.image.max()
import dask.multiprocessing, dask
RASTER_URL="https://github.com/mapbox/rasterio/raw/master/tests/data/rgb_deflate.tif"
items = [dask.delayed(worker)(RASTER_URL, roi) for roi in list_of_regions]
maxs = dask.compute(*items, get=dask.multiprocessing.get)
from numpy.ma.core import MaskedConstant
max([val for val in maxs if not isinstance(val,MaskedConstant)])
255