13

Let me state up front that I'm totally new to the field of GIS. I'm doing a project involving a database of buildings and their geometries in my home country (the Netherlands). At the moment, I'm trying to read geometry data from a geopackage using a python script.

If I understand the geopackage specification correctly, geometry data in a geopackage is a form of a well-known-binary. I tried to find a python package that can translate the binaries in the geopackage to usable geometries. I.e., I want to be able extract the points/linestrings/etc from the binaries.

I thought I could find a python package that would be able to do this, but I either don't understand how they work or they don't work with the data I have.

As an example, I have the following binary:

binary representation of geometry

or, as represented in python:

b'\x40\x71\x00\x00\x71\x3d\x0a\xd7\x05\xc3\x0e\x41\x23\xc3\x0e\x41\x71\x3d\x0a\xd7\xe3\x1c\x22\x41\xec\x1c\x22\x41\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\xeb\x03\x00\x00\x01\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00\x71\x3d\x0a\x41\xcd\xcc\xcc\xcc\xe5\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x18\xc3\x0e\x41\x71\x3d\x0a\x41\x00\x00\x00\x00\x00\x00\x00\x00\x29\x5c\x8f\x41\x52\xb8\x1e\x85\xea\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00'

When I try to read a binary using for example the shapely package, I get an error:

from shapely import wkb
geom = b'\x40\x71\x00\x00\x71\x3d\x0a\xd7\x05\xc3\x0e\x41\x23\xc3\x0e\x41\x71\x3d\x0a\xd7\xe3\x1c\x22\x41\xec\x1c\x22\x41\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\xeb\x03\x00\x00\x01\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00\x71\x3d\x0a\x41\xcd\xcc\xcc\xcc\xe5\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x18\xc3\x0e\x41\x71\x3d\x0a\x41\x00\x00\x00\x00\x00\x00\x00\x00\x29\x5c\x8f\x41\x52\xb8\x1e\x85\xea\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00'
point = wkb.loads(geom)
print(point)
shapely.errors.WKBReadingError: Could not create geometry because of errors while reading input.

To read geometry data from geopackages in python, will I have to write my own parser?

3
  • 2
    you are asking shapely to read a wkb geometry but your input is not in wkb format, this is why you are getting this error.
    – B-C B.
    Nov 24, 2019 at 22:21
  • Where did you get this binary from? This is not a valid GeoPackageBinaryHeader.
    – CL.
    Nov 25, 2019 at 14:30
  • @CL. It's this public dataset. Maybe I copied the sample wrong, but it should be a valid geopackage Nov 27, 2019 at 12:19

4 Answers 4

21

You can use the gdal/ogr, fiona (built on gdal/ogr) or geopandas (built on fiona) python libraries.

Below is a fiona example:

import fiona

# No need to pass "layer='etc'" if there's only one layer
with fiona.open('test.gpkg', layer='layer_of_interest') as layer:
    for feature in layer:
        print(feature['geometry'])

Partial output for one record in my data:

{'type': 'MultiPolygon', 'coordinates': [[[(147.01294051, -34.75046834699997), (147.01289292100003, -34.75075388199998), (147.01251220799998, -34.75068249799993), etc...

From the fiona docs:

A record you get from a collection is a Python dict structured exactly like a GeoJSON Feature. Fiona records are self-describing; the names of its fields are contained within the data structure and the values in the fields are typed properly for the type of record. Numeric field values are instances of type int and float, for example, not strings.

15

Geopandas is extremely useful and easy to use for this kind of thing: http://geopandas.org/

You can load vector data with attributes (e.g. from geopackages, shapefiles, etc.) into ?(geo)pandas dataframes which allow really easy analysis of the data without messing around with ogr/gdal.

import geopandas as gpd
data = gpd.read_file("path.mygeopackage.gpkg")
data.head()  # Prints the first 5 rows of the loaded data to see what it looks like.
0
9

GeoPackages are SQLite databases with a specific structure. You cannot just read geometries as "random" sequential bytes from a SQLite database, there might be fragmentation or similar.

If you want to do it low-level and without one of the fine libraries others suggested, the pure Python standard library way would be to use sqlite3 to open the file and then get the geometry bytes via SQL.

For example for http://www.geopackage.org/data/sample1_2.gpkg:

import sqlite3
conn = sqlite3.connect('sample1_2.gpkg')                                
cursor = conn.cursor()

# geometry column name is in gpkg_geometry_columns
# let'S pretend we already know it is 'Shape' for this dataset
cursor.execute("SELECT Shape FROM counties LIMIT 1;")
result = cursor.fetchone()[0]
# now you have the bytes of one geometry

It's not plain WKB though, you would still need to parse a binary header and then take the actual WKB bytes from the blob. I won't go into detail because it was obviously not what the questioner wanted, but if you want to do this the hard way, the specifications can be found at http://www.geopackage.org./spec/#gpb_format

4

if you have a geopackage file with extension .gpkg you can use ogr2ogr to translate it to a csv file with a geometry column like so:

ogr2ogr -f "CSV" MyLayer.csv MyLayer.gpkg -lco GEOMETRY=AS_WKT

If you use qgis you can do this manually using the convert format tool. You'll just need to input -lco GEOMETRY=AS_WKT in Additional creation options.

3
  • This already helps a bunch, thanks! Does this mean it's not (easily) possible to extract this from the geopackage file itself using a python library? I've already built the majority of a python script to extract all needed information from the geopackage. The geometries were the only field I was still having trouble with. Though I'm not unwilling, it would be a shame if I had to rewrite all that work. Nov 24, 2019 at 23:02
  • I am unsure I understand. Ogr2ogr uses Gdal python library. So I’d say it’s easy to extract this information from a gpkg using a python library. This code extracts all other information from the gpkg and extracts an extra column of geometry in wkt format. What is the output that you are looking for?
    – B-C B.
    Nov 25, 2019 at 2:49
  • Use the geopackage driver gdal.org/drivers/vector/gpkg.html from OGR in python or shell the command (subprocess.Popen(['ogr2ogr', '-f', '"esri shapefile"', os.path.join(os.environ.get('TEMP'),'MyLayer.shp'), 'MyLayer.gpkg', '-lco', 'GEOMETRY=AS_WKT'])) then open the shapefile to get the geometries. Nov 25, 2019 at 6:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.