项目需求
项目中需要使用影像做差值计算后根据其DN值进行密度分割渲染
项目构想
- 1、主要根据CreateColorRamp()添加条带
- 2、根据SetRasterColorTable()渲染
需要注意的是渲染值必须是int型或者byte,所以如果你的DN值是小数等,建议对数据进行拉伸处理,再进行渲染。
代码实现
- ReadTheRaster.py
def writeRasterInformation(self, data, Savepath, nband):
# driver = self.dataset.GetDriver()
driver = gdal.GetDriverByName('GTiff') # 自己指定生成类型
writeable = driver.Create(Savepath, self.cols, self.rows, self.bands, gdal.GDT_Byte) # TODO 新建数据集
writeable.SetGeoTransform(self.geotrans) # 写入仿射变换参数
writeable.SetProjection(self.proj) # 写入投影
for i in range(nband):
writeable.GetRasterBand(i + 1).WriteArray(data[i], 0, 0)
writeable.GetRasterBand(i + 1).SetNoDataValue(0) # todo 给各波段设置nodata值
writeable.GetRasterBand(i + 1).FlushCache() # todo 波段统计量
# 获取拉伸后数据最大值最小值 进行等分5分
mindata = writeable.GetRasterBand(i + 1).GetStatistics(0, 1)[0]
maxdata = writeable.GetRasterBand(i + 1).GetStatistics(0, 1)[1]
###渲染部分代码
itemdata = int((maxdata - mindata) / 5) # 等分五份取整
colortabledata = [(0, 255, 255), (72, 255, 72), (255, 255, 72), (255, 190, 72), (255, 72, 72), (255, 29, 29)]
ct = gdal.ColorTable()
for j in range(4):
ct.CreateColorRamp(int(mindata + j * itemdata), colortabledata[j], int(mindata + (j + 1) * itemdata),
colortabledata[j])
band = writeable.GetRasterBand(i + 1)
band.SetRasterColorTable(ct)
ct.CreateColorRamp(int(mindata + 4 * itemdata), colortabledata[4], int(maxdata), colortabledata[5])
band = writeable.GetRasterBand(i + 1)
band.SetRasterColorTable(ct)
###渲染代码
returnData=[mindata,maxdata,colortabledata]
print(writeable.GetRasterBand(i + 1).GetStatistics(0, 1)) # todo 计算波段统计量 输出为min\max \Mean\stddev
return returnData
- 调用
#TODO RASTERFUNCTION
def computoffset(ds,db): # TODO 计算偏移量从而得出行列号
data1=ds.computedoffset()
data2=db.computedoffset()
xoffset = int((data1[0] - data2[0]) / data1[1])
yoffset = int((data1[3] - data2[3]) / data1[5])
data = [xoffset, yoffset]
return data
def getcludedata(n,ds,db): # TODO 重新切片切成大小相同的
# data = db.getRasterInformation(n)[38:25697, 97:10678] #TODO 调整为多少行多少列 先化为相同行相同列 即二维数组切片
data = db.getRasterInformation(n)[computoffset(ds,db)[1]:ds.computeRows()+computoffset(ds,db)[1], computoffset(ds,db)[0]:db.computeCols()] # TODO 调整为多少行多少列 先化为相同行相同列 即二维数组切片
return data
def getcludedata2(n,ds,db):
data = ds.getRasterInformation(n)[0:ds.computeRows(),0:db.computeCols()-computoffset(ds,db)[0]]
return data
def getIntersection(n,ds,db,strainRate):
data1=getcludedata2(n,ds,db)
data2=getcludedata(n,ds,db)
data2[data1==0]=0 # todo 将data1中为0位置的元素对应把data2中的元素也置为0
data1[data2==0]=0 # todo 将data2中为0位置的元素对应把data1中的元素也置为0
result=abs(data2-data1)
return strainRate*result ###线性拉伸
ds = ReadTheRaster.ReadRaster(TheFirstfilepath) #TODO 文件路径
db = ReadTheRaster.ReadRaster(TheSecondefilepath)
status = ds.writeRasterInformation([getIntersection(1, ds, db,20)],outputfilepath, 1)