Python+GDAL练手(1)、读写栅格数据

import osgeo
from osgeo import ogr
from osgeo import gdal
from osgeo import osr
import os
import numpy as np
gdal.AllRegister()
osgeo.gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'NO')
osgeo.gdal.SetConfigOption('SHAPE_ENCODING', 'gb2312')

filePath = r'H:\task\mask_14.tif'#365个波段,25列,26行的tif
dataset = gdal.Open(filePath)
adfGeoTransform = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
nXSize = dataset.RasterXSize #列数
nYSize = dataset.RasterYSize #行数
im_bands = dataset.RasterCount #波段数
#print(im_bands)
data = dataset.ReadAsArray(0, 0, nXSize, nYSize)
#print(data)
##########以下操作属博主瞎搞,可略过
result=np.maximum(data,0)#判断小于0的像素值都赋值为0
print(result.shape)
test=result[10:80,0:24,0:25]
dataarray=np.mean(test,axis=0)对0轴方向取均值
print(dataarray.shape)
##########
driver = gdal.GetDriverByName("GTiff")
##一开始没有仔细看数据类型(gdal.GDT_Float32),导致出错,望下次注意
datasetnew = driver.Create( r'H:\task\new_80.tif', nXSize, nYSize, 1, gdal.GDT_Float32)
datasetnew.SetGeoTransform(adfGeoTransform)
datasetnew.SetProjection(im_proj)
band = datasetnew.GetRasterBand(1)
band.WriteArray(dataarray)









练习代码,勿喷,个人学习使用

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值