import os import glob from osgeo import gdal, gdalconst import subprocess import gdalTools import cv2 as cv import numpy as np from tqdm import tqdm from config_eval import ConfigEval import sys from Redundancy_predict import predict from Redundancy_predict_segmentation import predict_seg from osgeo import ogr,osr import os import gdalTools from shutil import rmtree, copyfile from skimage import morphology import argparse def polygonize(imagePath, raster_path, forest_shp_path): pwd = os.path.dirname(os.path.abspath(__file__)) os.chdir(os.path.realpath(os.path.join(pwd, 'polygonize/'))) polygonize_exe = os.path.realpath(os.path.join(pwd, 'polygonize/polygonize0529.exe')) # polygonize_exe = os.path.realpath(os.path.join(pwd, 'polygonize/polygonize1126.exe')) polygonize_path = os.path.realpath(os.path.join(pwd, 'polygonize/polygonize.config')) rmHole = "64" simpoly = "2" scale = "3" with open(polygonize_path,'w') as f_config: f_config.write("--image=" + imagePath+'\n') f_config.write("--edgebuf="+raster_path+'\n') f_config.write("--line="+forest_shp_path+'\n') f_config.write("--rmHole=" + rmHole + '\n') f_config.write("--simpoly=" + simpoly + '\n') f_config.write("--scale=" + scale) f_config.close() subprocess.call(polygonize_exe) def mkdir(path): if not os.path.exists(path): os.mkdir(path) ## python main.py --img D:\MyWorkSpace\myProject\SLP-CroplandExtraction\images\test.tif --weights ./ckpts/SOED3.pth --shp D:\MyWorkSpace\myProject\SLP-CroplandExtraction\results\result.shp if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--img', type=str, default=r'D:\2021\3\EESNet_test\test\test.tif', help='the path of image') parser.add_argument('--weights', type=str, default='D:/MyWorkSpace/paper/plough2/cropland_prediction/model_results/SOED3.pth', help='the path of weights') parser.add_argument('--shp', type=str, default=r'D:\MyWorkSpace\paper\plough\data\SOED2_cq\results.shp', help='the out path of shapefile') args = parser.parse_args() import time t1 = time.time() print('starting parsing arguments') cfg = ConfigEval() cfg.data_path = args.img cfg.model_path = args.weights outPath = args.shp outRoot = os.path.split(outPath)[0] outName = os.path.split(outPath)[-1].split('.')[0] print('starting inferring edge') predict(cfg) print('starting inferring binary map') predict_seg(cfg) out_shp_path = cfg.save_path + '_shp' image_root = cfg.data_path if os.path.exists(out_shp_path): rmtree(out_shp_path) mkdir(out_shp_path) mkdir(outRoot) imgList = [cfg.data_path] print('starting polygon') for imgPath in imgList: imageName = os.path.split(imgPath)[-1].split('.')[0] rasterPath = glob.glob(cfg.save_path + f'_ms/{imageName}*tif')[0] assert os.path.exists(imgPath), print(f'please cheack {imgPath}') assert os.path.exists(rasterPath), print(f'please cheack {rasterPath}') print(rasterPath) shpName = imageName + '_final.shp' shpPath = os.path.join(out_shp_path, shpName) out_skeleton_path = os.path.join(outRoot, 'skeleton.tif') im_proj, im_geotrans, im_width, im_height, edge_dn = gdalTools.read_img(rasterPath) ret2, th2 = cv.threshold(edge_dn, 0, 255, cv.THRESH_BINARY + cv.THRESH_OTSU) th2 = np.where(th2 > 0, 1, 0) skeleton = morphology.skeletonize(th2) skeleton = np.where(skeleton > 0, 255, 0).astype(np.uint8) gdalTools.write_img(out_skeleton_path, im_proj, im_geotrans, skeleton) polygonize(imgPath, out_skeleton_path, shpPath) assert os.path.exists(shpPath), f"please check your shapefile path:{shpPath}" rasterPath_seg = glob.glob(cfg.save_path + f'_seg/{imageName}*tif')[0] gdalTools.ZonalStatisticsAsTable(rasterPath_seg, shpPath) shp_baseName = os.path.basename(shpPath).split('.')[0] shp_root = os.path.split(shpPath)[0] shpList = glob.glob(f'{shp_root}/{shp_baseName}*') files = os.listdir(shp_root) for f in files: if f[-4:] == '.shp': shp_path = os.path.join(shp_root, f) ds = ogr.Open(shp_path, 0) layer = ds.GetLayer() layer.SetAttributeFilter("majority = 1 or majority = 2") driver = ogr.GetDriverByName('ESRI Shapefile') out_ds = driver.CreateDataSource(outPath) out_layer = out_ds.CopyLayer(layer, 'temp') del layer, ds, out_layer, out_ds polygonPath = outPath imgPath = rasterPath_seg linePath = os.path.join(outRoot, 'line.shp') out_polygon_path = os.path.join(outRoot, 'polygon.tif') out_line_path = os.path.join(outRoot, 'line.tif') out_line_dn_path = os.path.join(outRoot, 'line_dn.tif') gdalTools.pol2line(polygonPath, linePath) gdalTools.shp2Raster(linePath, imgPath, out_line_path, nodata=0) gdalTools.shp2Raster(polygonPath, imgPath, out_polygon_path, nodata=0) copyfile(rasterPath, out_line_dn_path) print(f'spend time:{time.time() - t1}s') os.system(f'python PatchRefinement.py')