[Py] 用 GDAL 將 Numpy Array 儲存為 GeoTiff

這邊參考網路上問答與 pyGDAL 官方網站的範例,將處理過的 Numpy Array 儲存為 GeoTiff,座標系統預設為 TWD 97 TM 121 ,可以自行更改(已寫入不同變數)。

 GeoTiff" ># -*- coding: utf-8 -*-
"""
Using pyGDAL to save Numpy Array to Raster
Created on Mon Nov 28 10:51:05 2016

@author: Sun Cheng-Wei
"""
from osgeo import gdal
import numpy as np
def Arr2Gtiff(NParray,x_size,y_size,x_min,y_max,dst_path):
    # TWD97 TM 121 = EPSG:3826    
    twd97_121_wkt_projection = '''PROJCS["TWD97 / TM2 zone 121",GEOGCS["TWD97",DATUM["Taiwan Datum 1997",SPHEROID["GRS 1980",6378137.0,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0.0,0.0,0.0,0.0,0.0,0.0,0.0],AUTHORITY["EPSG","1026"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","3824"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["central_meridian",121.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",250000.0],PARAMETER["false_northing",0.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3826"]]'''
    # TWD97 TM 119 = EPSG:3825
    twd97_119_wkt_projection = '''PROJCS["TWD97 / TM2 zone 119",GEOGCS["TWD97",DATUM["Taiwan Datum 1997",SPHEROID["GRS 1980",6378137.0,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0.0,0.0,0.0,0.0,0.0,0.0,0.0],AUTHORITY["EPSG","1026"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","3824"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["central_meridian",119.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",250000.0],PARAMETER["false_northing",0.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3825"]]'''
    # TWD67 121 TM 121 = EPSG:3828
    twd67_121_wkt_projection = '''PROJCS["TWD97 / TM2 zone 119",GEOGCS["TWD97",DATUM["Taiwan Datum 1997",SPHEROID["GRS 1980",6378137.0,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0.0,0.0,0.0,0.0,0.0,0.0,0.0],AUTHORITY["EPSG","1026"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","3824"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["central_meridian",119.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",250000.0],PARAMETER["false_northing",0.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3825"]]'''
    # TWD67 121 TM 119 = EPSG:3827
    twd67_119_wkt_projection = '''PROJCS["TWD67 / TM2 zone 119",GEOGCS["TWD67",DATUM["Taiwan Datum 1967",SPHEROID["GRS 1967 Modified",6378160.0,298.25,AUTHORITY["EPSG","7050"]],AUTHORITY["EPSG","1025"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","3821"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["central_meridian",119.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",250000.0],PARAMETER["false_northing",0.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3827"]]'''
    
    driver = gdal.GetDriverByName('GTiff')
    NParr_rev = NParray[::-1]
    outRaster = driver.Create(dst_path,NParr_rev.shape[1],NParr_rev.shape[0],1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((x_min, x_size, 0, y_max, 0, y_size))
    outRaster.SetProjection(twd97_121_wkt_projection)
    outBand = outRaster.GetRasterBand(1)
    outBand.WriteArray(NParr_rev)
    outBand.FlushCache()
    
    return outRaster
    
def Test():
    X_min = 248172
    Y_max = 2652130
    pixelWidth = 10
    pixelHeight = 10
    newRasterfn = './test.tif'
    array = np.array([[ 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, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 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, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
    Arr2Gtiff(array,pixelWidth,pixelHeight,X_min,Y_max,newRasterfn)

if __name__ == "__main__":
    Test()
GeoTiff 範例
GeoTiff 範例

請多多指教!

這個網站採用 Akismet 服務減少垃圾留言。進一步瞭解 Akismet 如何處理網站訪客的留言資料