可以把各个指数、光谱、以及火灾等级栅格组合成一个影像。把这个合成影像采样,
如何进行合成:进行波段组合,指数和光谱是来自于同一影像,这样行列就会保持一致,对于火灾等级栅格采样的时候要与对应的影像一致
网站下载数据与哥白尼等级对应不一致需要字段对应,与数据包里面的标准对上,才能当标签用也就是火灾等级得矢量形成栅格然后和就近时间得哨兵2图像对齐
处理完的结果是一个像素一个像素堆成
火灾等级矢量转为和输入栅格坐标、范围、行列一致的栅格,作为标签使用。
# ======================================================================
# 函数:矢量转栅格
# 参数:
# rasterInput:输入栅格(原始数据),作为 矢量转化成的栅格 的范围、投影等信息的参考
# shpFile:输入矢量数据
# fireClassField: 矢量中代表等级的字段名
# rasterOut:转化的栅格文件的输出路径,若为'',则不输出
#
# 输出:与输入数据行列、坐标一致的标签图
# ======================================================================
import numpy as np
from osgeo import gdal, gdalconst
from osgeo import ogr
import geopandas as gpd
def shp2raster(rasterInput, shpFile, fireClassField, rasterOut):
dataset = gdal.Open(rasterInput, gdalconst.GA_ReadOnly)
geo_transform = dataset.GetGeoTransform()
cols = dataset.RasterXSize # 列数
rows = dataset.RasterYSize # 行数
if isinstance(shpFile, str):
shp = ogr.Open(shpFile,0)
print('\n输入矢量为shp的路径字符串')
else:
print(type(shpFile.to_json()))
shp = ogr.Open(shpFile.to_json())
print('\n输入矢量为GeoDataFrame格式的SHP')
m_layer = shp.GetLayer()
if rasterOut: # 判断是否进行输出,若输出则默认为tif格式,且rasterOut为完整路径
target_ds = gdal.GetDriverByName('GTiff').Create(rasterOut, cols,
rows, 1,
gdal.GDT_Byte)
print("tif")
else: # 存入内存,不进行输出,且rasterOut为''
target_ds = gdal.GetDriverByName('MEM').Create(rasterOut, cols,
rows,1,
gdal.GDT_Byte)
print("MEM")
target_ds.SetGeoTransform(geo_transform)
print('geo_transform:',geo_transform)
band = target_ds.GetRasterBand(1)
# band.SetNoDataValue(0) # 不执行该句后,背景值变成0;执行后背景值为Nodata (nan)
band.FlushCache()
# # gdal.RasterizeLayer(ds, bands, layer, burn_values, options=["BURN_VALUE_FROM=Z"])
gdal.RasterizeLayer(target_ds, [1], m_layer, options=["ATTRIBUTE="+fireClassField]) # 跟shp字段给栅格像元赋值
# crs = osr.SpatialReference()
# crs.ImportFromEPSG(raster_epsg_number)
# target_ds.SetProjection(crs.ExportToWkt()) # !!!!
target_ds.SetProjection(dataset.GetProjection()) # python3.6版本,投影要在矢量化后,推测是因为矢量的范围比栅格的范围要小,先投影会出错(3.9不会有问题)
shp.Release()
# return dataset, target_ds
if __name__ == '__main__':
rasterInput = 'E:\Raster\EMSR132_post.tif' # 原影像
rasterOut = "E:\Raster\EMSR132_lable_test.tif"
shpFile = gpd.read_file("E:\SHP\shp_EMSR132_fire_values.shp") # 原始输入矢量shp
tc_shp_area = shpFile.area
# print(tc_shp_area)
print('原有的总矢量特征数:', len(shpFile))
for i in range(len(tc_shp_area)):
if np.isnan(tc_shp_area[i]) or (tc_shp_area[i] == 0):
shpFile.drop([i], inplace=True)
print('去掉空白矢量特征后的总特征数:', len(shpFile))
bbb = shp2raster(rasterInput, shpFile,"grading" ,rasterOut)
del bbb
把标签合并到特征影像上,你可以用代码写,就是标签数组和特征影像的数组的合并,然后保存成新的栅格。
采样的话,如果是用机器学习,例如随机森林(一般用sk-learning库)等,需要把栅格采样成一个个点,作为分类器的输入;
如果是深度学习的话,需要裁剪把栅格裁剪成一个个影像块,作为其输入。
代码实现
RGB三通道、外加近红、短波红外、NBR、NDVI、BAI、dNBR11个通道需要自己拿波段出来的数据进行实验
多找几个火灾,切成小块,尽可能多点好 做一下敏感性分析,找到适合分类的最理想的就是云量、烟雾、亮度
遇到数据集不丰富的问题
考虑做进行波段组合,就是说可以同一个shp文件,随机波段合并,然后拿到指数的敏感度,也可以作为训练集,(但是我不知道不敏感的数据会不会影响效果,但是一般大家用的都是公认敏感的数据),然后标签数组和特征影像的数组进行合并,就有新的栅格数据,然后把栅格采样成一个个影像块,就扩大了数据集,然后输入
如需yolo模型,需要tif转图像,不然不能输入训练
但是tif文件提出来的三通道jpg格式全黑,rgb提出来然后裁成小格,但是那个火灾矢量文件提rgb出来全黑
统计一下那个全黑图像,看看是否全是0,里面一般会有有1、2、3、4,分别代表不同等级,是否等级全是0?
最后失败
可能需要专业遥感文件
关于ARCGIS推荐码农设计师公众号 内容丰富,全面
最后还是失败了,处理的图像不能跑yolo