个人整理的python地理信息库,本文代码也在其中，欢迎使用: https://github.com/Adam0429/geo-py
顺便帮公司打个广告，欢迎遥感类专家的加入～～
投简历戳这里，英视睿达欢迎您！
废话不多说，直接上代码
from tqdm import tqdm
from osgeo import gdal,ogr,osr
import numpy as np
from glob import glob
def read_tif(tif_path):
ds = gdal.Open(tiff_path)
row = ds.RasterXSize
col = ds.RasterYSize
band = ds.RasterCount
for i in range(band):
data = ds.GetRasterBand(i+1).ReadAsArray()
data = np.expand_dims(data , 2)
if i == 0:
allarrays = data
else:
allarrays = np.concatenate((allarrays, data), axis=2)
return {'data':allarrays,'transform':ds.GetGeoTransform(),'projection':ds.GetProjection(),'bands':band,'width':row,'height':col}
# 左上角点坐标 GeoTransform[0],GeoTransform[3] Transform[1] is the pixel width, and Transform[5] is the pixel height
def write_tif(fn_out, im_data, transform,proj=None):
'''
功能:
----------
将矩阵按某种投影写入tif，需指定仿射变换矩阵，可选渲染为rgba
参数:
----------
fn_out:str
输出tif图的绝对文件路径
im_data: np.array
tif图对应的矩阵
transform: list/tuple
gdal-like仿射变换矩阵，若im_data矩阵起始点为左上角且投影为4326，则为
(lon_x.min(), delta_x, 0,
lat_y.max(), 0, delta_y)
proj: str（wkt格式）
投影，默认投影坐标为4326，可用osr包将epsg转化为wkt格式，如
srs = osr.SpatialReference()# establish encoding
srs.ImportFromEPSG(4326) # WGS84 lat/lon
proj = srs.ExportToWkt() # create wkt fromat of proj
# 设置投影，proj为wkt format
if proj is None:
proj = 'GEOGCS["WGS 84",\
DATUM["WGS_1984",\
SPHEROID["WGS 84",6378137,298.257223563, \
AUTHORITY["EPSG","7030"]], \
AUTHORITY["EPSG","6326"]], \
PRIMEM["Greenwich",0, \
AUTHORITY["EPSG","8901"]], \
UNIT["degree",0.0174532925199433, \
AUTHORITY["EPSG","9122"]],\
AUTHORITY["EPSG","4326"]]'
# 渲染为rgba矩阵
# 设置数据类型
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 将(通道数、高、宽)顺序调整为(高、宽、通道数)
# print('shape of im data:', im_data.shape)
im_bands = min(im_data.shape)
im_shape = list(im_data.shape)
im_shape.remove(im_bands)
im_height, im_width = im_shape
band_idx = im_data.shape.index(im_bands)
# 找出波段是在第几个
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(fn_out, im_width, im_height, im_bands, datatype)
# if dataset is not None:
dataset.SetGeoTransform(transform) # 写入仿射变换参数
dataset.SetProjection(proj) # 写入投影
if im_bands == 1:
# print(im_data[:, 0,:].shape)
if band_idx == 0:
dataset.GetRasterBand(1).WriteArray(im_data[0, :, :])
elif band_idx == 1:
dataset.GetRasterBand(1).WriteArray(im_data[:, 0, :])
elif band_idx == 2:
dataset.GetRasterBand(1).WriteArray(im_data[:, :, 0])
else:
for i in range(im_bands):
if band_idx == 0:
dataset.GetRasterBand(i + 1).WriteArray(im_data[i, :, :])
elif band_idx == 1:
dataset.GetRasterBand(i + 1).WriteArray(im_data[:, i, :])
elif band_idx == 2:
dataset.GetRasterBand(i + 1).WriteArray(im_data[:, :, i])
dataset.FlushCache()
del dataset
driver = None
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
使用案例
>>> geo.read_tif('/Users/wangfeihong/Desktop/GF2_PMS2_E102.4_N26.3_20210420_L1A0005608668-MSS2_clip.tif')
{'data': array([[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
...,
[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
...,
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]]], dtype=uint16), 'transform': (102.36514511007022, 7.796940285642904e-06, 0.0, 26.32195006415771, 0.0, -7.797041881219687e-06), 'projection': 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]', 'bands': 4, 'width': 11577, 'height': 6836}
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50