如何在python中編寫/創建GeoTIFF RGB圖像文件？

Question

我有形狀爲nx 5個numpy的陣列，NY

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

的紅色，綠色和藍色陣列包含從0-255範圍和緯度/經度陣列是緯度值/經度像素座標。

我的問題是如何將這些數據寫入geotiff？

我最終想用底圖繪製圖像。

這是我到目前爲止的代碼，但是我得到一個巨大的GeoTIFF文件（〜500MB），它出現空白（只是一個黑色的圖像）。還要注意的是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

我認爲這個問題是當你創建的數據集，你通過它GDT_Float32。對於像素範圍爲0-255的標準圖像，您需要GDT_Byte。 Float要求值通常在0-1之間。

我把你的代碼和隨機生成的一些數據，所以我可以測試你的API的其餘部分。然後，我在太浩湖周圍創建了一些虛擬座標。

這是代碼。

#!/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 

這裏是輸出。 （注：我們的目標是產生色彩，不受地形！）

這裏是QGIS的圖像，驗證投影。