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.
I don't have a lot of aerial imagery to test this on just now, so I'm not 100% sure it'll work in your case. But it should get you closer to what you want.
I've used the rasterio library, which provides a very nice wrapper around GDAL and makes the code much simpler to write (and read).
This combines several things:-
- a generator function which passes back the name of each TIFF file in a folder (and any sub-folders).
- it opens each band in each file to generate a list of booleans, which represents which bands have NODATA values
- if all bands have NODATA values, write to the console
Code :-
import os
import rasterio
def find_images(path):
for root, dirs, files in os.walk(path, topdown=False):
for name in files:
if os.path.splitext(name)[-1] == ".tiff":
yield os.path.join(root, name)
def main():
for file_name in find_images("/path/to/your/images"): # <--- change this!
with rasterio.open(file_name, "r", driver="GTiff") as source:
no_data = source.nodata
bands = source.read()
got_nulls = [no_data in band for band in bands]
if False not in got_nulls: # we want list elements to all be True
print("File {} NodataValue {}, has {} Bands, Got Nulls {}".format(file_name, no_data, len(bands), got_nulls))
if __name__ == "__main__":
main()
I've tried this on Python 2.7.12 using rasterio 0.26 and numpy 1.8.2.
As with most python experiments, I recommend using a python virtual environment to avoid having to install anything system-wide :)
Best Answer
1- One way is contacting the developer/owner/project manager/scientist who puts the
.tif
files online.2- Another way is to read the
.tif
description (some kind of internal metadata that is available via.tif
format.3- Is to locate the bounds and look at its aerial/satellite images and speculate what it might be.
4- If your raster is accompanied by an
*.hdr
orreadme.txt
, open them in a text viewer and you'll possibly find all the information