如何在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的图像，验证投影。