Import raster into PostGIS using Python

Since this doesn't have an answer, and I don't quite have an appreciation for the PostGIS documentation (yet?), I'll post my solution.

For me, I needed to get GeoTIFFs in programmatically. Calling an external program wasn't preferable (for project relevant reasons), but seems to be the default for PostGIS interaction.

PostGIS has support for a few GDAL compatible raster types, and can translate between them freely (well at the expense of cheap CPU cycles, at least). The hard part is getting something GDAL compatible into PsycoPG2 (especially without GDAL). If you can make a GeoTIFF of the data you're interested in (using OpenCV or something else) then you've got a decent amount of the puzzle solved.

def load_into_PostGIS(connection, image_name):
    with open(image_name, 'rb') as f:
        with connection: # To autocommit/rollback
            with connection.cursor() as cursor:
                cursor.execute("INSERT INTO table(rast) VALUES (ST_FromGDALRaster(%s))", (f.read(),))

PostGIS should be able to get all the raster information from the GeoTIFF or other compatible image format just fine. If you have an image in memory, not on disk, it'll probably be easier to save it to disk, but you can create an in memory file-like bytestring using OpenCV at least with cv2.imencode or with the GDAL library with a little more work. I'm a little fuzzy on how to get a numpy image into GDAL since that isn't in my workflow, but if you need help, I can look into it. (It should be easier than saving the image to disk, opening it with GDAL in update mode, then adding GeoReference data, but that would be a workaround if you need it.) As for getting the bytestring from GDAL, it ended up being as easy as:

# Solution found buried in a mailing list exchange:
# https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html
def bytes_from_gdal(ds):
    # Don't worry, this makes it in RAM, not on disk.
    # We need the virtual file support to be able to read it, afaik.
    vsipath = '/vsimem/output.tif'
    gdal.GetDriverByName("GTiff").CreateCopy(vsipath, ds)

    # Read /vsimem/output.tif
    f = gdal.VSIFOpenL(vsipath, 'rb')
    gdal.VSIFSeekL(f, 0, 2) # seek to end
    size = gdal.VSIFTellL(f)
    gdal.VSIFSeekL(f, 0, 0) # seek to beginning
    data = gdal.VSIFReadL(1, size, f)
    gdal.VSIFCloseL(f)
    return data

Which you can directly put into the raster column:

def insert_into_postgis(connection, byte):
    with connection: # To autocommit/rollback
        with connection.cursor() as cursor:
            cursor.execute("INSERT INTO table(rast) VALUES (ST_FromGDALRaster(%s))", (f.read(),))

So if you have a BAG file loaded in gdal as ds and a psycopg2 connection opened as conn, you just call the above as:

insert_into_postgis(conn, bytes_from_gdal(ds))

And for sickening completeness, the opposite direction is much easier. You can use the answer here for that. Once you get the bytes from the database, you can save the image to disk as a binary file and load it as a GeoTIFF however you normally would. Changing it to BAG is another issue, as I believe BAG is a readonly file format for GDAL currently (correct me if you can).

There's also likely a more direct way to do it, rather than using GeoTIFFs as an intermediate, but I'm not a GIS guy. Someone more familiar with this software should help out if they can. I know I would appreciate it.