Comment écrire/créer un fichier image GeoTIFF RVB en python?

Question

I ont 5 matrices de numpy de nx ny, de forme

lons.shape = (nx,ny) 
lats.shape = (nx,ny) 
reds.shape = (nx,ny) 
greens.shape = (nx,ny) 
blues.shape = (nx,ny) 

Les rouges, verts et bleus des tableaux de valeurs contiennent allant de 0 à 255 et les matrices de latitude/longitude sont pixel coordonnées latitude/longitude.

Ma question est comment puis-je écrire ces données à un géotiff? Je souhaite finalement tracer l'image en utilisant le fond de carte.

Voici le code que j'ai jusqu'à présent, mais je reçois un énorme fichier GeoTIFF (~ 500 Mo) et il apparaît vide (juste une image noire). Notez également que nx, ny = 8120, 5416.

from osgeo import gdal 
from osgeo import osr 
import numpy as np 
import h5py 
import os 

os.environ['GDAL_DATA'] = "/Users/andyprata/Library/Enthought/Canopy_64bit/User/share/gdal" 

# read in data 
input_path = '/Users/andyprata/Desktop/modisRGB/' 
with h5py.File(input_path+'red.h5', "r") as f: 
    red = f['red'].value 
    lon = f['lons'].value 
    lat = f['lats'].value 

with h5py.File(input_path+'green.h5', "r") as f: 
    green = f['green'].value 

with h5py.File(input_path+'blue.h5', "r") as f: 
    blue = f['blue'].value 

# convert rgbs to uint8 
r = red.astype('uint8') 
g = green.astype('uint8') 
b = blue.astype('uint8') 

# set geotransform 
nx = red.shape[0] 
ny = red.shape[1] 
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()] 
xres = (xmax - xmin)/float(nx) 
yres = (ymax - ymin)/float(ny) 
geotransform = (xmin, xres, 0, ymax, 0, -yres) 

# create the 3-band raster file 
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 3, gdal.GDT_Float32) 
dst_ds.SetGeoTransform(geotransform) # specify coords 
srs = osr.SpatialReference()   # establish encoding 
srs.ImportFromEPSG(3857)    # WGS84 lat/long 
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file 
dst_ds.GetRasterBand(1).WriteArray(r) # write r-band to the raster 
dst_ds.GetRasterBand(2).WriteArray(g) # write g-band to the raster 
dst_ds.GetRasterBand(3).WriteArray(b) # write b-band to the raster 
dst_ds.FlushCache()      # write to disk 
dst_ds = None       # save, close 

Answer 1

Je pense que le problème est lorsque vous créez le Dataset, vous passez GDT_Float32. Pour les images standard avec des plages de pixels de 0 à 255, vous avez besoin de GDT_Byte. Float nécessite des valeurs entre 0 et 1 généralement.

J'ai pris votre code et généré aléatoirement des données afin que je puisse tester le reste de votre API. J'ai ensuite créé quelques coordonnées factices autour du lac Tahoe.

Voici le code.

#!/usr/bin/env python 
from osgeo import gdal 
from osgeo import osr 
import numpy as np 
import os, sys 

# Initialize the Image Size 
image_size = (400,400) 

# Choose some Geographic Transform (Around Lake Tahoe) 
lat = [39,38.5] 
lon = [-120,-119.5] 

# Create Each Channel 
r_pixels = np.zeros((image_size), dtype=np.uint8) 
g_pixels = np.zeros((image_size), dtype=np.uint8) 
b_pixels = np.zeros((image_size), dtype=np.uint8) 

# Set the Pixel Data (Create some boxes) 
for x in range(0,image_size[0]): 
    for y in range(0,image_size[1]): 
     if x < image_size[0]/2 and y < image_size[1]/2: 
      r_pixels[y,x] = 255 
     elif x >= image_size[0]/2 and y < image_size[1]/2: 
      g_pixels[y,x] = 255 
     elif x < image_size[0]/2 and y >= image_size[1]/2: 
      b_pixels[y,x] = 255 
     else: 
      r_pixels[y,x] = 255 
      g_pixels[y,x] = 255 
      b_pixels[y,x] = 255 

# set geotransform 
nx = image_size[0] 
ny = image_size[1] 
xmin, ymin, xmax, ymax = [min(lon), min(lat), max(lon), max(lat)] 
xres = (xmax - xmin)/float(nx) 
yres = (ymax - ymin)/float(ny) 
geotransform = (xmin, xres, 0, ymax, 0, -yres) 

# create the 3-band raster file 
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 3, gdal.GDT_Byte) 

dst_ds.SetGeoTransform(geotransform) # specify coords 
srs = osr.SpatialReference()   # establish encoding 
srs.ImportFromEPSG(3857)    # WGS84 lat/long 
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file 
dst_ds.GetRasterBand(1).WriteArray(r_pixels) # write r-band to the raster 
dst_ds.GetRasterBand(2).WriteArray(g_pixels) # write g-band to the raster 
dst_ds.GetRasterBand(3).WriteArray(b_pixels) # write b-band to the raster 
dst_ds.FlushCache()      # write to disk 
dst_ds = None 

Voici la sortie. (Note: Le but est de produire des couleurs, pas térrain!)

Voici l'image dans QGIS, la validation de la projection.