南通网站建设空间专业建设网站

接上一部分继续介绍~

一、地类矢量转栅格

这一步是为了能让地类值和影像的对象落在同一区域，从而将影像中的分割对象同化为实际地物类别。

train_fn = r".\train_data1.shp"
train_ds = ogr.Open(train_fn)
lyr = train_ds.GetLayer()
driver = gdal.GetDriverByName('MEM')
target_ds = driver.Create('', im_width, im_height, 1, gdal.GDT_UInt16)
target_ds.SetGeoTransform(im_geotrans)
target_ds.SetProjection(im_proj)
options = ['ATTRIBUTE=tyPE']
gdal.RasterizeLayer(target_ds, [1], lyr, options=options)
data = target_ds.GetRasterBand(1).ReadAsArray()
ground_truth = target_ds.GetRasterBand(1).ReadAsArray()
ground_truth = ground_truth.transpose((1, 0))
classes = np.unique(ground_truth)[1:]

最终得到带有地物类别数据的栅格点数据。

二、特征匹配

将得到的栅格点真实地物数据通过迭代与影像对象相匹配后，通过迭代器寻找对象的相应特征。

segments_per_class = {}
for klass in classes:segments_of_class = segments[ground_truth == klass]segments_per_class[klass] = set(segments_of_class)intersection = set()
accum = set()
for class_segments in segments_per_class.values():intersection |= accum.intersection(class_segments)accum |= class_segments
assert len(intersection) == 0, "Segment(s) represent multiple classes"
print("adjust complete")
end3 = time.time()
print('Running time2: %s Seconds'%(end3-start3))print("start train randomforest classification")
start4 = time.time()
train_img = np.copy(segments)
threshold = train_img.max() + 1 for klass in classes:class_label = threshold + klassfor segment_id in segments_per_class[klass]:train_img[train_img == segment_id] = class_labeltrain_img[train_img <= threshold] = 0
train_img[train_img > threshold] -= thresholdtraining_objects = []
training_labels = []
for klass in classes:class_train_object = [v for i, v in enumerate(objects) if segment_ids[i] in segments_per_class[klass]]training_labels += [klass] * len(class_train_object)training_objects += class_train_object

在实际的影像样本构建过程中，有的地物样本可能彼此距离相差较小，造成两个或多个样本落在同一个分割区域，这样会导致特征匹配迭代无限进行下去，所以我们要从两个或多个样本中取其一。

三、分类

最后利用scikit-learn的随机森林分类器，对样本分割块和其他未定义分割块进行预测，最后将结果输出到栅格中。

def PolygonizeTheRaster(inputfile,outputfile):dataset = gdal.Open(inputfile, gdal.GA_ReadOnly)srcband=dataset.GetRasterBand(1)im_proj = dataset.GetProjection()prj = osr.SpatialReference() prj.ImportFromWkt(im_proj)drv = ogr.GetDriverByName('ESRI Shapefile')dst_ds = drv.CreateDataSource(outputfile)dst_layername = 'out'dst_layer = dst_ds.CreateLayer(dst_layername, srs=prj)dst_fieldname = 'DN'fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)dst_layer.CreateField(fd)dst_field = 0gdal.Polygonize(srcband, None, dst_layer, dst_field) classifier = RandomForestClassifier(n_jobs=-1)  
classifier.fit(training_objects, training_labels) 
predicted = classifier.predict(objects)  clf = segments.copy()
for segment_id, klass in zip(segment_ids, predicted):clf[clf == segment_id] = klass
# temp = temp.transpose((2, 1, 0))
mask = np.sum(temp, axis=2)  
mask[mask > 0.0] = 1.0   
mask[mask == 0.0] = -1.0clf = np.multiply(clf, mask)
clf[clf < 0] = -9999.0
clf = clf.transpose((1, 0))
clfds = driverTiff.Create(r"D:\Data\testdata\classification\result.tif", im_width, im_height,1, gdal.GDT_Float32) 
clfds.SetGeoTransform(im_geotrans)
clfds.SetProjection(im_proj)
clfds.GetRasterBand(1).SetNoDataValue(-9999.0)
clfds.GetRasterBand(1).WriteArray(clf)
clfds = Noneend4 = time.time()
outputfile = r".\result2.shp"
print('Running time2: %s Seconds'%(end4-start4))
PolygonizeTheRaster(r".\result.tif",outputfile)

类在影像中的标记ID是1、2、3等等数值，因此用一般的栅格查看软件很难从肉眼进行查看。这里为了方便读者查看以及制图，我还进行了栅格转矢量的操作，这样放到arcmap中能清晰的查看分类情况。

最后的分类结果：

图1 分类结果图

图2 林地分类效果