How To Overwrite A GeoPackage Raster

This came to me as somewhat of a shock, I was bounty hunting at gis.stackexchange
when I came across this question and took a look at the GDAL formats docs.
Apparently, you can't overwrite a raster stored inside a GeoPackage. It's simply not one of the options. But you can however overwrite an existing vector layer with since it's one of the creation options.

So I decided to test it out, if it was really impossible, or just not implemented.
First thing I did was simply look what tables were inside a GeoPackage, and since I'm too lazy to read the standard and the schema I just opened one with good old DBeaver.

Apparently, there are quite a lot of tables in a GeoPackage...

 

But, did you notice?
A lot of those tables are from the three vector layers inside, so I just looked everywhere else.
And then I started deleting every mention of the `DEM` raster layer (made a point to write down everywhere it was).

 

What I came up with was the following code, which enables you to erase and then replace a raster layer from inside a GeoPackage, If you know the name of the layer.


    import os
    import sqlite3
    gpkg = 'path/to/your/db.gpkg'
    layer_name= 'superawesomeraster'


    conn = sqlite3.connect(gpkg)
    conn.execute('DROP TABLE IF EXISTS {}'.format(layer_name))
    conn.execute('DELETE FROM gpkg_2d_gridded_coverage_ancillary WHERE tile_matrix_set_name = "{}"'.format(layer_name))
    conn.execute('DELETE FROM gpkg_2d_gridded_tile_ancillary WHERE tpudt_name="{}"'.format(layer_name));
    conn.execute('DELETE FROM gpkg_tile_matrix WHERE table_name="{}"'.format(layer_name));
    conn.execute('DELETE FROM gpkg_tile_matrix_set WHERE table_name="{}"'.format(layer_name));
    conn.execute('DELETE FROM gpkg_contents WHERE table_name="{}"'.format(layer_name));
    conn.execute('DELETE FROM gpkg_extensions WHERE table_name="{}"'.format(layer_name));
    conn.commit()
    conn.execute('VACUUM')

    conn.close()
 
    inputRas = layer_in.source()
    outputRas = gpkg
    gdal_string = 'gdal_translate -of GPKG "{}" "{}" -co RASTER_TABLE={} -co APPEND_SUBDATASET=YES'.format(inputRas, outputRas,lyr) 
    os.system(gdal_string) 
 

What this does, is use a simple connection to the GeoPackage (which can be connected to like any other spatialite db) and delete all mention of the table with basic SQL.
Then simply use gdal_translate to write another layer in, in this case a QGIS layer. This part can be replaced with any other way of writing rasters into a GeoPackage, but I find gdal_translate the simplest.

Comments