J'ai des données raster dans une table postgres que je veux entrer dans python en tant que tableau numpy. J'utilise psycopg2 pour me connecter à la base de données. Je peux télécharger les données mais elles reviennent sous forme de chaîne (probablement un binaire sérialisé).

Est-ce que quelqu'un sait comment prendre cette chaîne et la convertir en un tableau numpy?

J'ai exploré d'autres options pour télécharger le raster comme utiliser st_astiff et encoder pour télécharger le fichier hex et utiliser xxd mais cela n'a pas fonctionné. Je continue à obtenir l'erreur «rt_raster_to_gdal: Impossible de charger le pilote GDAL de sortie» et je n'ai pas les autorisations pour définir les variables d'environnement pour pouvoir activer les pilotes.

TL, DR: souhaitez importer des données raster dans un tableau numpy (en utilisant python).

rt_raster_to_gdal: impossible de charger le pilote GDAL de sortie

Quant à la première erreur avec ST_AsTIFF , vous devez activer vos pilotes GDAL, qui par défaut ne sont pas activés pour PostGIS 2.1. Consultez le manuel pour savoir comment procéder. Par exemple, j'ai une variable d'environnement configurée sur un ordinateur Windows avec:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

qui peut être confirmé avec PostGIS avec:

SELECT short_name, long_name
FROM ST_GDALDrivers();

PostGIS à Numpy

Vous pouvez exporter la sortie vers un fichier GeoTIFF de mémoire virtuelle pour que GDAL puisse le lire dans un tableau Numpy. Pour obtenir des conseils sur les fichiers virtuels utilisés dans GDAL, consultez cet article de blog .

import os
import psycopg2
from osgeo import gdal

# Adjust this to connect to a PostGIS database
conn = psycopg2.connect(...)
curs = conn.cursor()

# Make a dummy table with raster data
curs.execute("""\
    SELECT ST_AsRaster(ST_Buffer(ST_Point(1, 5), 10), 10, 10, '8BUI', 1) AS rast
    INTO TEMP mytable;
""")

# Use a virtual memory file, which is named like this
vsipath = '/vsimem/from_postgis'

# Download raster data into Python as GeoTIFF, and make a virtual file for GDAL
curs.execute("SELECT ST_AsGDALRaster(rast, 'GTiff') FROM mytable;")
gdal.FileFromMemBuffer(vsipath, bytes(curs.fetchone()[0]))

# Read first band of raster with GDAL
ds = gdal.Open(vsipath)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()

# Close and clean up virtual memory file
ds = band = None
gdal.Unlink(vsipath)

print(arr)  # this is a 2D numpy array

Affiche un point tamponné tramé.

[[0 0 0 1 1 1 1 0 0 0]
 [0 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 0]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [0 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 0]
 [0 0 0 1 1 1 1 0 0 0]]

Notez que j'ai utilisé un format «GTiff» dans l'exemple, mais d' autres formats pourraient être mieux adaptés. Par exemple, si vous avez un grand raster qui doit être transféré via une connexion Internet lente, essayez d'utiliser 'PNG' pour le compresser.

Je pense que la question était de savoir si vous pouvez lire à partir de tables raster postgis SANS pilotes gdal activés. Comme tout ce qui concerne Python, vous le pouvez!

Assurez-vous de sélectionner votre résultat raster comme WKBinary:

sélectionnez St_AsBinary (rast) ...

Utilisez le script ci-dessous pour déchiffrer WKBinary dans un format d'image python. Je préfère l'opencv, car il gère un nombre arbitraire de bandes d'images, mais on peut utiliser PIL / low si 1 ou 3 bandes sont plus habituelles.

Je ne gère que les images d'octets pour l'instant, mais il est relativement trivial de les étendre à d'autres types de données.

J'espère que c'est utile.

structure d'importation
import numpy as np
importer cv2

# Fonction pour déchiffrer l'en-tête WKB
def wkbHeader (brut):
    # Voir http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

    en-tête = {}

    header ['endianess'] = struct.unpack ('B', raw [0]) [0]
    header ['version'] = struct.unpack ('H', raw [1: 3]) [0]
    header ['nbands'] = struct.unpack ('H', raw [3: 5]) [0]
    header ['scaleX'] = struct.unpack ('d', raw [5:13]) [0]
    header ['scaleY'] = struct.unpack ('d', raw [13:21]) [0]
    header ['ipX'] = struct.unpack ('d', raw [21:29]) [0]
    header ['ipY'] = struct.unpack ('d', raw [29:37]) [0]
    header ['skewX'] = struct.unpack ('d', raw [37:45]) [0]
    header ['skewY'] = struct.unpack ('d', raw [45:53]) [0]
    header ['srid'] = struct.unpack ('i', raw [53:57]) [0]
    header ['width'] = struct.unpack ('H', raw [57:59]) [0]
    header ['height'] = struct.unpack ('H', raw [59:61]) [0]

    en-tête de retour

# Fonction pour déchiffrer les données raster WKB 
def wkbImage (brut):
    h = wkbHeader (brut)
    img = [] # array pour stocker les bandes d'images
    offset = 61 # longueur brute d'en-tête en octets
    pour i dans la plage (h ['nbands']):
        # Déterminez le type de pix pour ce groupe
        pixtype = struct.unpack ('B', raw [offset]) [0] >> 4
        # Pour l'instant, nous ne gérons que les octets non signés
        si pixtype == 4:
            band = np.frombuffer (raw, dtype = 'uint8', count = h ['width'] * h ['height'], offset = offset + 1)
            img.append ((np.reshape (band, ((h ['height'], h ['width'])))))
            offset = offset + 2 + h ['largeur'] * h ['hauteur']
        # à faire: gérer d'autres types de données 

    return cv2.merge (tuple (img))