问题描述
有M个study area, 每一个study area有多种类别,每一种类别有N个polygons,要求从每个study area的每一类中采样得到N个点。
解决办法
- 采样可以用arcMap
- 批量处理可以用python遍历文件夹
相关
- kml to shapefile: ogrmerge.py
- setup python and ArcGIS
- 可以直接在arcmap的python窗口里面load代码然后运行
- 如果用pycharm运行python,需要设置project interpreter
Sample code
1 convert kml files to shapefiles
#!/bin/bash
set -e
#the dir of the kml files
Dir='/data/qiu/LCZ_wudpat_train/raw/*/'
files=$(find $Dir -name '*.km*')
for f in $files; do
echo $f
city0=${f%/*}
#cut f: delete the first (starting from right) / and everything rightwards.
# echo $city0
city=${city0##*/}
#cut city0: delete the first (starting from left) / and everything leftwards.
echo $city
rm -r -f $city0/$city'_shapeFile'
rm -r -f /data/qiu/LCZ_wudpat_train/shapeFile/$city
ogrmerge.py -o /data/qiu/LCZ_wudpat_train/shapeFile/$city $f -nln $city'_'{LAYER_NAME}
#the output can be one file or multiple files for multiple subdir in the kml file
done;
# Ref:
# https://gdal.org/programs/ogrmerge.html
2 sampling with arcpy
import os
import arcpy
import numpy as np
import shutil
import glob
#top_folder = r"F:\0Data4Paper\LCZ_wudpat_train\test"
#'folder that contains multiple folders for study areas'
top_folder = r"F:\0Data4Paper\LCZ_wudpat_train\shapeFile"
#'tmp dir for intermediate results'
intermFolder = r"F:\0Data4Paper\LCZ_wudpat_train\tmp"
for path, dirs, files in os.walk(top_folder):
for d in dirs:
#print(path, dirs, files)
sh_path = os.path.join(top_folder, d)
#print(sh_path)
#'''clear existing files'''
saveDisFolder=os.path.join(r'F:\0Data4Paper\LCZ_wudpat_train\ds', d)
print(saveDisFolder)
if os.path.exists(saveDisFolder):
shutil.rmtree(saveDisFolder)
if not os.path.exists(saveDisFolder):
os.mkdir(saveDisFolder)
pointsFolder = os.path.join(r"F:\0Data4Paper\LCZ_wudpat_train\points", d)
print(pointsFolder)
if os.path.exists(pointsFolder):
shutil.rmtree(pointsFolder)
if not os.path.exists(pointsFolder):
os.mkdir(pointsFolder)
if os.path.exists(intermFolder):
shutil.rmtree(intermFolder)
if not os.path.exists(intermFolder):
os.mkdir(intermFolder)
#for each class:
for i in np.arange(1,18):
c=i
if i>10:
i=i+90
file=os.path.join(sh_path, d+'_'+str(i) + "." + 'shp')
print(file)
if os.path.isfile(file):
#dissolve polygons into one and save temporally
saveDis=os.path.join(saveDisFolder, d+'_'+str(i))
arcpy.Dissolve_management(file, saveDis, "", "", "", "")
#sample points from the dissolved polygons
arcpy.CreateRandomPoints_management(intermFolder, d+'_'+str(i), saveDis+ "." + 'shp', "", 30, "100 Meter", "MULTIPOINT")
savePoints=os.path.join(intermFolder, d+'_'+str(i))
#add a field to assign classes
arcpy.AddField_management(savePoints+ "." + 'shp', "class", "SHORT")
arcpy.CalculateField_management(savePoints+ "." + 'shp', "class", int(c))
else:
print('maybe this class not exist!')
fileList = glob.glob(intermFolder+ "\*.shp")
# #print(fileList)
#'merge points into one file for each study area'
arcpy.Merge_management(fileList, os.path.join(pointsFolder, d), "")