將PilotGaea double matrix 轉成 geotifff

array2tiff.py

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 14 09:03:09 2020

@author: User
"""

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

# config
GDAL_DATA_TYPE = gdal.GDT_Float32
GEOTIFF_DRIVER_NAME = r'GTiff'
NO_DATA = -999999
SPATIAL_REFERENCE_SYSTEM_WKID = 4326  

def create_raster(output_path,
                  columns,
                  rows,
                  nband = 1,
                  gdal_data_type = GDAL_DATA_TYPE,
                  driver = GEOTIFF_DRIVER_NAME):
    ''' returns gdal data source raster object

    '''
    # create driver
    driver = gdal.GetDriverByName(driver)

    output_raster = driver.Create(output_path,
                                  int(columns),
                                  int(rows),
                                  nband,
                                  eType = gdal_data_type)    
    return output_raster

def numpy_array_to_raster(output_path,
                          numpy_array,
                          upper_left_tuple,
                          cell_resolution,
                          nband = 1,
                          no_data = NO_DATA,
                          gdal_data_type = GDAL_DATA_TYPE,
                          spatial_reference_system_wkid = SPATIAL_REFERENCE_SYSTEM_WKID,
                          driver = GEOTIFF_DRIVER_NAME):
    ''' returns a gdal raster data source

    keyword arguments:

    output_path -- full path to the raster to be written to disk
    numpy_array -- numpy array containing data to write to raster
    upper_left_tuple -- the upper left point of the numpy array (should be a tuple structured as (x, y))
    cell_resolution -- the cell resolution of the output raster
    nband -- the band to write to in the output raster
    no_data -- value in numpy array that should be treated as no data
    gdal_data_type -- gdal data type of raster (see gdal documentation for list of values)
    spatial_reference_system_wkid -- well known id (wkid) of the spatial reference of the data
    driver -- string value of the gdal driver to use

    '''

    print('UL: (%s, %s)' % (upper_left_tuple[0],upper_left_tuple[1]))

    rows, columns = numpy_array.shape
    print('ROWS: %s\n COLUMNS: %s\n' % (rows,columns))

    # create output raster
    output_raster = create_raster(output_path,
                                  int(columns),
                                  int(rows),
                                  nband,
                                  gdal_data_type) 

    
    geotransform = (upper_left_tuple[0],
                    cell_resolution,
                    0,
                    upper_left_tuple[1] + cell_resolution,
                    0,
                    cell_resolution*-1)    
    
    

    spatial_reference = osr.SpatialReference()
    spatial_reference.ImportFromEPSG(spatial_reference_system_wkid)
    output_raster.SetProjection(spatial_reference.ExportToWkt())
    output_raster.SetGeoTransform(geotransform)
    output_band = output_raster.GetRasterBand(1)
    output_band.SetNoDataValue(no_data)
    output_band.WriteArray(numpy_array)          
    output_band.WriteArray(numpy_array)
    output_band.WriteArray(numpy_array)
    output_band.FlushCache()
    output_band.ComputeStatistics(False)

    if os.path.exists(output_path) == False:
        raise Exception('Failed to create raster: %s' % output_path)

    return  output_raster


local_vars = {}

def test():
    global local_vars
    output_path = "D:\\dem.tif"
    
    dem = open("D:\\dem.txt", "r")
    
    
    line = dem.readline()
    line = dem.readline()
    
    Rows=748
    Columns=748
    
    numpy_array  = np.empty([Rows,Columns], dtype=float)
    cnt = 0    
    line = dem.readline()
    while(line):
        lines = line.split(" ")
        row = np.asarray(lines, dtype=np.float32)
        numpy_array[Rows-cnt-1,] = row
        cnt += 1
        line = dem.readline()
    
    upper_left_tuple = np.array([119.324956,24.481719])
    cell_resolution = (123.067800-119.324956)/748
    numpy_array_to_raster(output_path,numpy_array,upper_left_tuple,cell_resolution)
    local_vars = inspect.currentframe().f_locals
    
    output_path = "D:\\vm.tif"
    
    dem = open("D:\\vm.txt", "r")
    
    
    line = dem.readline()
    line = dem.readline()
    
    Rows=748
    Columns=748
    
    numpy_array  = np.empty([Rows,Columns], dtype=float)
    cnt = 0    
    line = dem.readline()
    while(line):
        lines = line.split(" ")
        row = np.asarray(lines, dtype=np.float32)
        numpy_array[Rows-cnt-1,] = row
        cnt += 1
        line = dem.readline()
    
    upper_left_tuple = np.array([119.324956,24.481719])
    cell_resolution = (123.067800-119.324956)/748
    numpy_array_to_raster(output_path,numpy_array,upper_left_tuple,cell_resolution)
    local_vars = inspect.currentframe().f_locals    
    
test()  

source : https://gis.stackexchange.com/questions/290776/how-to-create-a-tiff-file-using-gdal-from-a-numpy-array-and-specifying-nodata-va?noredirect=1&lq=1

precondition : conda install -c conda-forge gdal