[GIS] How to read MrSID files in Python (using gdal)

gdalmrsidpython

I have to read a large number of MrSID images and process them in Python. I usually access GIS file formats using the gdal application and the gdal Python library. However, the gdal package for Python does not come with MrSID support (I use the portable WinPython and the precompiled gdal package from here). Then I also use QGIS that comes with gdal support: here, gdal does support MrSID out of the box (command: gdalinfo.exe –formats).

I was wondering if there is a gdal package for Python that does support MrSID out of the box, or do I have to compile it from the source (using Windows 8.1 64-bit).

I tried to just copy the QGIS dll "…\gdalplugins\gdal_MrSID.dll" to the same folder in my Python gdal library, but this makes the command "gdalinfo.exe –formats" to stop with an error.

This is the code that I use in Python to see if gdal supports MrSID:

import gdal

gdal.AllRegister()
for i in range(1, gdal.GetDriverCount()):
    drv = gdal.GetDriver(i)
    print drv.GetDescription(),

Using this after implementing Kersten's solution, the list does not include the entry "MrSID". Calling the command "gdalinfo –formats" however, does list the entry "MrSID (rov): Multi-resolution Seamless Image Database (MrSID)".

Best Answer

Instead of using the binaries provided by Christopher Gohlke you can use the GDAL binaries from GISInternals:

For your case you would need:

Related Solutions

[GIS] 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.

[GIS] Building GDAL with libcurl support on Debian Linux

Both of those commands are wrong. It should be something like,

./configure --with-curl=/usr/bin/curl-config

curl-config is an executable which provides details on the curl version and locations of libraries, etc.

You can get this information by typing,

./configure --help|grep curl

But as EvilGenius says, you shouldn't need this command at all. Make sure curl is actually installed, of course.

MB Edit:

The correct libcurl installation syntax is:

apt-get install libcurl4-gnutls-dev

(once I had that installed, the config option above worked great and the installation allowed the CartoDB format to show up)

Best Answer

Related Solutions

[GIS] Import raster into PostGIS using Python

[GIS] Building GDAL with libcurl support on Debian Linux

Related Question