本文幾乎全部內容來自博客 https://www.cnblogs.com/ninicwang/p/11533066.html
1、gdal包簡介
gdal是空間數據處理的開源包，其支持超過100種柵格數據類型，涵蓋所有主流GIS與RS數據格式。
2、讀取柵格數據
#導入gdal包
from osgeo import gdal
#導入numpy包（支持高維數組和矩陣運算，也提供了許多數組和矩陣運算的函數）
import numpy as np
#打開文件
dataset=gdal.Open("fdem.tif")
#柵格矩陣的列數（X是列）
im_width = dataset.RasterXSize
#柵格矩陣的行數（Y是行）
im_height = dataset.RasterYSize
#波段數
im_bands = dataset.RasterCount
#共有六個參數，分表代表左上角x坐標；東西方向上圖像的分辨率；如果北邊朝上，地圖的旋轉角度，0表示圖像的行與x軸平行；左上角y坐標；
im_geotrans = dataset.GetGeoTransform()
>>> im_geotrans
(409294.88696681266, 27.376482012944024, 0.0, 4423871.083377095, 0.0, -27.376482012944006)
#地圖投影信息
im_proj = dataset.GetProjection()
>>> im_proj
'PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1
984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG
","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AU
THORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUT
HORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION[
"Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PAR
AMETER["central_meridian",117],PARAMETER["scale_factor",0.99
96],PARAMETER["false_easting",500000],PARAMETER["false_north
ing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easti
ng",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]'
#讀取某一像素點的值
#（1）讀取一個波段，其參數為波段的索引號，波段索引號從1開始(我打開的這幅圖像只有一個波段)
band=dataset.GetRasterBand(1)
#（2）用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，讀出從(xoff,yoff)開始，大小為(xsize,ysize)的矩陣。以下為讀取整幅圖像
im_datas=band.ReadAsArray(0,0,im_width,im_height)
#（3）獲取某一或某幾個像素的值(查看10~14 行和 20~25 列的數據)
data=im_datas[10:15,20:26]
#釋放內存。如果不釋放，在arcgis或envi中打開該圖像時顯示文件已被占用
del dataset