My guess is that the gdal-dev version you have is outdated and that the version of python binding you got in pip is not compatible any more. I accidentally uninstalled python-gdal from RHEL, and unsuccessfully tried to install gdal with pip today. The first error I got was
extensions/gdal_wrap.cpp:3373: error: ‘CPLGetErrorHandlerUserData’ was not declared in this scope
extensions/gdal_wrap.cpp: In function ‘CPLErr PushErrorHandler(void (*)(CPLErr, int, const char*), void*)’:
which is the same you got near the beginning. I downloaded gdal v1.10.1 source today, and I grepped CPLGetErrorHanlerUserData
, and found it in include/cpl_error.h, and not in the file in version 1.8.1 I had. After I compiled gdal 1.10.1, installed it and rerun pip install gdal
. It gave a bit of warning but installed successfully.
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.
Best Answer
If you want GDAL you can use this (I found this code in another post and it works nice):
I prefer use rasterio with fiona to get the coordinates and use the sample tool:
Both solutions are just for points!