读取tif

'''

# from libtiff import TIFF

import matplotlib.pyplot as plt

#导入numpy包（支持高维数组和矩阵运算，也提供了许多数组和矩阵运算的函数）

from osgeo.gdal import GDT_Float32

import numpy as np

from osgeo import gdal

# path =

#打开文件

dataset=gdal.Open(r"/data/publicData/hxfLUCC/GLC_FCS30_2020_E100N30.tif")

#栅格矩阵的列数

im_width = dataset.RasterXSize

#栅格矩阵的行数

im_height = dataset.RasterYSize

#波段数

im_bands = dataset.RasterCount

#仿射矩阵，左上角像素的大地坐标和像素分辨率。

#共有六个参数，分表代表左上角x坐标；东西方向上图像的分辨率；如果北边朝上，地图的旋转角度，0表示图像的行与x轴平行；左上角y坐标；

#如果北边朝上，地图的旋转角度，0表示图像的列与y轴平行；南北方向上地图的分辨率。

im_geotrans = dataset.GetGeoTransform()

print(im_geotrans)

#地图投影信息

im_proj = dataset.GetProjection()

print(im_proj)

#读取某一像素点的值

#（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]

print(im_datas)

plt.imshow(im_datas)

plt.show()

#释放内存。如果不释放，在arcgis或envi中打开该图像时显示文件已被占用

del dataset

driver = gdal.GetDriverByName("GTiff")

#创建空文件，并确定开辟多大内存；每个像素都有一个对应的值，这个值得类型用数据类型指定。这里的数据类型是gdal数据类型。

dataset = driver.Create("./hjnTEST.tif", im_datas.shape[1],im_datas.shape[0], 1, GDT_Float32)

#设置头文件信息

#（1）写入仿射变换参数

dataset.SetGeoTransform(im_geotrans)

#（2）写入投影信息(这里我所用的投影是从上一篇文章里的tif文件读来的)

dataset.SetProjection(im_proj)

#写入数据体

#为了方便，这里的im_bands是从上一篇文章里的tif文件所读取到的波段数，im_data是从上一篇文章里的tif文件所读取到的数据

im_bands=1

if im_bands == 1:

    dataset.GetRasterBand(1).WriteArray(im_datas) #写入数组数据

else:

    for i in range(im_bands):

        dataset.GetRasterBand(i+1).WriteArray(im_datas[i])

#释放内存空间

del dataset

'''