赞
踩
废话不多说,直接上代码
读取tif:
- from osgeo import gdal
- import sys
- def Read_img2array(img_file_path):
- """
- 读取栅格数据,将其转换成对应数组
- img_file_path: 栅格数据路径
- :return: 返回投影,几何信息,和转换后的数组
- """
- dataset = gdal.Open(img_file_path) # 读取栅格数据
- print('处理图像波段数总共有:', dataset.RasterCount)
- # 判断是否读取到数据
- if dataset is None:
- print('Unable to open *.tif')
- sys.exit(1) # 退出
- projection = dataset.GetProjection() # 投影
- geotrans = dataset.GetGeoTransform() # 几何信息
- im_width = dataset.RasterXSize #栅格矩阵的列数
- im_height = dataset.RasterYSize #栅格矩阵的行数
- im_bands = dataset.RasterCount #波段数
- # 直接读取dataset
- img_array = dataset.ReadAsArray()
- return im_width,im_height,im_bands,projection, geotrans, img_array
数组写入tif:
- import numpy as np
- from osgeo import gdal
- #tiff_file为tif,im_data为数组
- def Write_img2array(tiff_file, im_proj, im_geotrans, data_array):
- if 'int8' in data_array.dtype.name:
- datatype = gdal.GDT_Int16
- elif 'int16' in data_array.dtype.name:
- datatype = gdal.GDT_Int16
- else:
- datatype = gdal.GDT_Float32
-
- if len(data_array.shape) == 3:
- im_bands, im_height, im_width = data_array.shape
- else:
- im_bands, (im_height, im_width) = 1,data_array.shape
-
- driver = gdal.GetDriverByName("GTiff")
- dataset = driver.Create(tiff_file, im_width, im_height, im_bands, datatype)
-
- dataset.SetGeoTransform(im_geotrans)
- dataset.SetProjection(im_proj)
-
- if im_bands == 1:
- dataset.GetRasterBand(1).WriteArray(data_array)
- else:
- for i in range(im_bands):
- dataset.GetRasterBand(i+1).WriteArray(data_array[i])
-
- del dataset
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。