How to query stored PostGIS Raster data from within the IPython Notebook

Question

Test data

• 1 TIF file (159KB)

Goal threefold:

1. Load raster into PostGIS using raster2pgsql and visualize in QGIS
2. In my IPython Notebook connect to PostGIS and load raster into NumPy array
3. In my IPyhton Notebook use Pandas to load a time-series of one pixel of rasters with different time step stored in PostGIS

Process so far I've managed to get one raster image into PostGIS Raster using the raster2pgsql command and visualize it in QGIS using the DB Manager:

raster2pgsql -s 4326 -d -I -C -M -R -l 4 D:\Downloads\raster//modis.tif -F -t 100x100 public.ndvi | psql -U postgres -d rastertest -h localhost -p 5432

But how to access/query this raster from within IPython Notebook?

I found this presentation, which is going about SQLALchemy and GeoAlchemy2. And where it is mentioned that it support PostGIS Raster as well. It seems to be very interesting! But using the documentation I don't see how I can apply this to Raster data

I think I can make a connection to my PostGIS database using the following code, where postgres=user, password=admin and database=rastertest:

from sqlalchemy import create_engine
engine = create_engine('postgresql://postgres:admin@localhost/rastertest', echo=True)

But then.. any advice is very much appreciated.

Answer 1

You should use the psycopg module to connect to a postgres db from python. Some code samples:

import psycopg2


def connect_db():

try:
    conn = psycopg2.connect("dbname='rastertest' user='admin' host='localhost' password='password'")
    conn.set_session(autocommit=True) #if you want your updates to take effect without being in a transaction and requiring a commit, for a beginner, I would set this to True
    return conn
except:
    print "I am unable to connect to the database"
    return None

def get_raster(raster_id,conn):

    query= "SELECT ST_AsText(geom) from raster_table where id={}".format(raster_id)
    conn.cursor().execute(query)
    res = cur.fetchall()

    return res[0][0]

Maybe the text representation of the raster is something you can use. Alternatively, take a look here http://postgis.net/docs/RT_reference.html to see if any of the functions return what you want for your numpy array and replace the query in get_raster accordingly. (Possibly this http://postgis.net/docs/RT_ST_DumpValues.html)