命令行记录-python gdal读取栅格数据

本文几乎全部内容来自博客 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