问题描述

有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), "")