基于对象的土地覆盖分类：将 shapefile 转换为栅格

如何解决基于对象的土地覆盖分类：将 shapefile 转换为栅格

对于我的硕士论文，我目前正在尝试创建一个基于对象的分类代码。在创建我的分割形状文件并光栅化我的训练数据后，我陷入困境并搜索如何在代码中实现我的 shapefile。我想尝试栅格化 shapefile，但它包含统计信息（17 个属性），因此属性应该以某种方式变成带。以下是我的代码中棘手的部分目前的样子：

# Filename of the raster Tiff that will be created
raster_fn = f'{classif_path}object.tif'

# Filename of input OGR file
train_fn = f'{shapefile_path}SegmentationSTAT_BCK_median.shp'
train_ds = ogr.Open(train_fn)
lyr = train_ds.GetLayer()
# create a new raster layer in memory
driver = gdal.GetDriverByName('MEM')
target_ds = driver.Create('',median_temp_tif.RasterXSize,median_temp_tif.RasterYSize,1,gdal.GDT_UInt16)
target_ds.SetGeoTransform(median_temp_tif.GetGeoTransform())
target_ds.SetProjection(median_temp_tif.GetProjection())
# rasterize the training points
options = ['ATTRIBUTE=A2017_7__m','ATTRIBUTE=A2017_7__1','ATTRIBUTE=A2017_7__2','ATTRIBUTE=A2017_7__s','ATTRIBUTE=A2018_7__m','ATTRIBUTE=A2018_7__1','ATTRIBUTE=A2018_7__2','ATTRIBUTE=A2018_7__s','ATTRIBUTE=A2019_7__m','ATTRIBUTE=A2019_7__1','ATTRIBUTE=A2019_7__2','ATTRIBUTE=A2019_7__s','ATTRIBUTE=A2020_7__m','ATTRIBUTE=A2020_7__1','ATTRIBUTE=A2020_7__2','ATTRIBUTE=A2020_7__s','ATTRIBUTE=CODE']
gdal.RasterizeLayer(target_ds,[17],lyr,options=options)
# retrieve the rasterized data and print basic stats
data = target_ds.GetRasterBand(1).ReadAsArray()

这里是整个分类代码：

in_situ_cal_shp = f'{shapefile_path}{period}_cal.shp' 
in_situ_cal_tif = f'{classif_path}{period}_cal.tif' 
median_temp_tif = f'{median_path}S1A_IW_SLC__1SDV_A2018_7_DESCENDING_Orb_Thermal_Cal_Sigma_Deb_TC_SRTM1s_Clip_AF_VH_filtered_mean_norm_medianFilter9x9.tif'
polygones_val_shp = f'{shapefile_path}{period}_val.shp'
polygones_val_tif = f'{classif_path}{period}_val.tif'


# Open the calibration shapefile with GeoPandas

in_situ_gdf = gpd.read_file(in_situ_cal_shp)

# Open the raster file you want to use as a template for rasterize

img_src = rasterio.open(median_temp_tif)
img     = img_src.read()

meta = img_src.meta.copy()
meta.update(compress='lzw')
meta.update(nodata=0)

# Burn the features into the raster and write it out

dst_src = rasterio.open(in_situ_cal_tif,'w+',**meta)
dst_arr = dst_src.read(1)

# this is where we create a generator of geom,value pairs to use in rasterizing

geom_col = in_situ_gdf.geometry
code_col = in_situ_gdf.CODE

shapes = ((geom,value) for geom,value in zip(geom_col,code_col))


in_situ_cal_arr = features.rasterize(shapes=shapes,fill=0,out=dst_arr,transform=dst_src.transform)

dst_src.write_band(1,in_situ_cal_arr)

# Close rasterio objects

img_src.close()
dst_src.close()

# Filename of the raster Tiff that will be created
raster_fn = f'{classif_path}object.tif'

# Filename of input OGR file
train_fn = f'{shapefile_path}SegmentationSTAT_BCK_median.shp'
train_ds = ogr.Open(train_fn)
lyr = train_ds.GetLayer()
# create a new raster layer in memory
driver = gdal.GetDriverByName('MEM')
target_ds = driver.Create('',options=options)
# retrieve the rasterized data and print basic stats
data = target_ds.GetRasterBand(1).ReadAsArray()

# Create an empty list to append all feature rasters one by one
list_src_arr = []

list_im = target_ds

for im_file in list_im:

    src = rasterio.open(im_file,"r")
    im = src.read(1)
    list_src_arr.append(im)
    src.close()

# Merge all the 2D matrices from the list into one 3D matrix

feat_arr = np.dstack(list_src_arr)

print(feat_arr.shape)
print(f'There are {feat_arr.shape[2]} features')
print(f'The features type is : {feat_arr.dtype}')

# Open in-situ used for calibration

src = rasterio.open(in_situ_cal_tif,"r")
cal_arr = src.read(1)
src.close()


# Find how many non-zero entries we have -- i.e. how many training data samples?

n_samples = (cal_arr > 0).sum()

print(f'We have {n_samples} samples\n')


# What are our classification labels?
labels = np.unique(cal_arr[cal_arr > 0])
print(f'The training data include {labels.size} classes: {labels} \n')


# We will need a "X" matrix containing our features,and a "y" array containing our labels
#     These will have n_samples rows
#     In other languages we would need to allocate these and them loop to fill them,but NumPy can be faster

X = feat_arr[cal_arr > 0,:]
y = in_situ_cal_arr[cal_arr > 0]

X = np.nan_to_num(X,nan=-10)

print(f'Our X matrix is sized: {X.shape}')
print(f'Our y array is sized: {y.shape}')

start_training = time.time()

# Initialize our model with 500 trees
rf = RandomForestClassifier(n_estimators=500,oob_score=True)

# Fit our model to training data
rf = rf.fit(X,y)

end_training = time.time()

# Get time elapsed during the Random Forest training
hours,rem = divmod(end_training-start_training,3600)
minutes,seconds = divmod(rem,60)
print("Random Forest training : {:0>2}:{:0>2}:{:05.2f}".format(int(hours),int(minutes),seconds))

print(f'Our OOB prediction of accuracy is: {round(rf.oob_score_ * 100,2)}%')

bands = range(1,X.shape[1]+1)

for b,imp in zip(bands,rf.feature_importances_):
    print(f'Band {b} importance: {round(imp,4)}')

# Setup a dataframe
df = pd.DataFrame()

df['truth'] = y
df['predict'] = rf.predict(X)

# Cross-tabulate predictions

cross_tab = pd.crosstab(df['truth'],df['predict'],margins=True)

display(cross_tab)

# Take our full image and reshape into long 2d array (nrow * ncol,nband) for classification

img = feat_arr

img = np.nan_to_num(img,nan=-10)

new_shape = (img.shape[0] * img.shape[1],img.shape[2])

img_as_array = img[:,:,:].reshape(new_shape)

print(f'Reshaped from {img.shape} to {img_as_array.shape}')


start_classification = time.time()

# Now predict for each pixel
class_prediction = rf.predict(img_as_array)

# Reshape our classification map
class_prediction = class_prediction.reshape(img[:,0].shape)

end_classification = time.time()

hours,rem = divmod(end_classification-start_classification,seconds))

classif_tif = f'{classif_path}Classif_{period}.tif'

print(class_prediction.shape)

img_src = rasterio.open(median_temp_tif)
img = img_src.read()

meta = img_src.meta.copy()
meta.update(compress='lzw')

print(meta)


dst = rasterio.open(classif_tif,'w',**meta)

dst.write(class_prediction,1)

dst.close()




# Open the shapefile with GeoPandas

polygones_gdf = gpd.read_file(polygones_val_shp)

# Open the raster file you want to use as a template for rasterize

img_src = rasterio.open(median_temp_tif)
img     = img_src.read()

meta = img_src.meta.copy()
meta.update(compress='lzw')
meta.update(nodata=0)

# Burn the features into the raster and write it out

dst_src = rasterio.open(polygones_val_tif,value pairs to use in rasterizing

geom_col = polygones_gdf.geometry
code_col = polygones_gdf.CODE

shapes = ((geom,code_col))


polygones_val_arr = features.rasterize(shapes=shapes,polygones_val_arr)

# Close rasterio objects

img_src.close()
dst_src.close()