API tools faq
paste
Advertisement
Standa

CHM_heights_tops

Standa
Aug 3rd, 2015
257
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 15.83 KB | None | 0 0
raw download clone embed print report
  1. #1 import modulu z ArcGIS
  2. import arcpy
  3.  
  4. from os import listdir
  5. from os.path import isfile, join
  6.  
  7. #2 kontrola extenzi
  8. arcpy.CheckOutExtension("3D")
  9. arcpy.CheckOutExtension("spatial")
  10.  
  11. arcpy.env.overwriteOutput = True
  12. arcpy.env.outputCoordinateSystem = 102067 # to by mel byt kod pro S-JTSK_Krovak_East_North
  13.  
  14. #3 class Preprocessing: vyber vstupnich souboru a nastaveni ulozist i pro ostatni tridy,
  15. # orez dat podle hranic SLP, prevod .las do .shp (tzv. Multipointu)
  16. class PreProcessing():
  17. def __init__(self):
  18. self.first = "98_any_first.las", "99_any_first.las", "100_any_first.las", "101_any_first.las"
  19. self.first_dir = "C:\\Data LLS pro IGA\\data_IGA_klasif\\any_first_classified\\"
  20. self.first_res = "C:\\cviceni_Python\\uloziste_MTP\\MTP_af_3b.shp"
  21. self.first_res_clip = "C:\\cviceni_Python\\uloziste_MTP\\MTP_af_3b_clip.shp"
  22. self.ground = "98_gr_mod.las", "99_gr_mod.las", "100_gr_mod.las", "101_gr_mod.las"
  23. self.ground_dir = "C:\\Data LLS pro IGA\\data_IGA_klasif\\ground_default_classified\\"
  24. self.ground_res = "C:\\cviceni_Python\\uloziste_MTP\\MTP_gr_3b.shp"
  25. self.ground_res_clip = "C:\\cviceni_Python\\uloziste_MTP\\MTP_gr_3b_clip.shp"
  26. self.SLP_hranice = "C:\\Data LLS pro IGA\\LLS_zajmove\\Hranice_Vrstevnice_SLP\\slp_hranice.shp"
  27. self.cestaPlochy = "C:\\cviceni_Python\\plochy_2014_fce_3b\\"
  28. self.vystupniAdresar = "C:\\cviceni_Python\\uloziste\\"
  29. self.f_dir_drop = "C:\\cviceni_Python\\uloziste\\flow_dir_drop"
  30. self.VystupniAdresarReclassify = "C:\\Users\\Stanislav\\Documents\\ArcGIS\\Default4.gdb\\"
  31. self.VystupniAdresarVYSLEDKY = "C:\\cviceni_Python\\uloziste_vysledky\\"
  32. self.VystupniAdresarEXCEL = "C:\\cviceni_Python\\uloziste_vysledky_excel\\"
  33.  
  34. def zdrojove_MTP(self, a, s):
  35. return [a+x for x in s]
  36.  
  37. def do_preprocessing(self):
  38. print(self.zdrojove_MTP(self.first_dir, self.first))
  39. arcpy.LASToMultipoint_3d(\
  40. self.zdrojove_MTP(self.first_dir, self.first), \
  41. self.first_res, "0.5", "", "ANY_RETURNS", "", \
  42. "PROJCS['S-JTSK_Krovak_East_North',GEOGCS['GCS_S_JTSK',DATUM['D_S_JTSK',SPHEROID['Bessel_1841',6377397.155,299.1528128]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Krovak'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Pseudo_Standard_Parallel_1',78.5],PARAMETER['Scale_Factor',0.9999],PARAMETER['Azimuth',30.28813975277778],PARAMETER['Longitude_Of_Center',24.83333333333333],PARAMETER['Latitude_Of_Center',49.5],PARAMETER['X_Scale',-1.0],PARAMETER['Y_Scale',1.0],PARAMETER['XY_Plane_Rotation',90.0],UNIT['Meter',1.0]]", "las", "1", "NO_RECURSION")
  43.  
  44. print(self.zdrojove_MTP(self.ground_dir, self.ground))
  45. arcpy.LASToMultipoint_3d(\
  46. self.zdrojove_MTP(self.ground_dir, self.ground), \
  47. self.ground_res, "0.5", "", "ANY_RETURNS", "", \
  48. "PROJCS['S-JTSK_Krovak_East_North',GEOGCS['GCS_S_JTSK',DATUM['D_S_JTSK',SPHEROID['Bessel_1841',6377397.155,299.1528128]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Krovak'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Pseudo_Standard_Parallel_1',78.5],PARAMETER['Scale_Factor',0.9999],PARAMETER['Azimuth',30.28813975277778],PARAMETER['Longitude_Of_Center',24.83333333333333],PARAMETER['Latitude_Of_Center',49.5],PARAMETER['X_Scale',-1.0],PARAMETER['Y_Scale',1.0],PARAMETER['XY_Plane_Rotation',90.0],UNIT['Meter',1.0]]", "las", "1", "NO_RECURSION")
  49.  
  50. arcpy.Clip_analysis(self.first_res, self.SLP_hranice, self.first_res_clip, "")
  51. arcpy.Clip_analysis(self.ground_res, self.SLP_hranice, self.ground_res_clip, "")
  52.  
  53. def nastavPlochy(self):
  54. vstupniSoubory=[]
  55. for f in listdir(self.cestaPlochy):
  56. if isfile(join(self.cestaPlochy, f)):
  57. if (f[-4:]==".shp"):
  58. vstupniSoubory.append(join(self.cestaPlochy, f))
  59. return vstupniSoubory
  60.  
  61. def nastavVystup(self, poradi, prvni="_pl_af"):
  62. vystupniSoubory=[]
  63. for i in range(1, poradi + 1):
  64. vystupniSoubory.append(self.vystupniAdresar+str(i)+prvni)
  65. return vystupniSoubory
  66.  
  67. #4 class CHMPreparation: priprava Canopy height modelu, orez dat pouze na kruhove zajmove plochy,
  68. # tvorba DMT a DMP, jejich odecteni a primarni zhlazeni nizkofrekvencni filtraci pro dalsi pouziti
  69. class CHMPreparation():
  70. def __init__(self):
  71. self.nastavPlochy = pp.nastavPlochy()
  72. self.first_res_clip = pp.first_res_clip
  73. self.ground_res_clip = pp.ground_res_clip
  74. self.VystupPrvni = pp.nastavVystup(len(pp.nastavPlochy()), "_pl_af.shp")
  75. self.VystupPosledni = pp.nastavVystup(len(pp.nastavPlochy()), "_pl_gr.shp")
  76. self.VystupIDW_af = pp.nastavVystup(len(pp.nastavPlochy()), "_idw_af")
  77. self.VystupIDW_gr = pp.nastavVystup(len(pp.nastavPlochy()), "_idw_gr")
  78. self.VystupMINUS = pp.nastavVystup(len(pp.nastavPlochy()), "_minus")
  79. self.VystupFILTER_LOW = pp.nastavVystup(len(pp.nastavPlochy()), "_filter_low")
  80.  
  81. def do_CHMPreparation(self):
  82. for i in range(0, len(self.nastavPlochy)):
  83. arcpy.Clip_analysis(self.first_res_clip, self.nastavPlochy[i], self.VystupPrvni[i])
  84. arcpy.Clip_analysis(self.ground_res_clip, self.nastavPlochy[i], self.VystupPosledni[i])
  85.  
  86. arcpy.Idw_3d(self.VystupPrvni[i], "Shape.Z", self.VystupIDW_af[i], "0.5", "2", "VARIABLE 12", "")
  87. arcpy.Idw_3d(self.VystupPosledni[i], "Shape.Z", self.VystupIDW_gr[i], "0.5", "2", "VARIABLE 12", "")
  88.  
  89. arcpy.Minus_3d(self.VystupIDW_af[i], self.VystupIDW_gr[i], self.VystupMINUS[i])
  90.  
  91. arcpy.gp.Filter_sa(self.VystupMINUS[i], self.VystupFILTER_LOW[i], "LOW", "DATA")
  92.  
  93. #5 class CHMProcessing: filtrace CHM 8-mi variantami nastaveni Focal Statistics
  94. # (zasadni vliv na finalni pocet vrcholu, mensi vliv na mereni vysek)
  95. # vynasobeni vsech variant "-1" pro dalsi zpracovani CHM pomoci nastroju hydrologickeho modelovani
  96. class CHMProcessing():
  97. def __init__(self):
  98. self.VystupFILTER_LOW = chmpre.VystupFILTER_LOW
  99. self.vystupniAdresar = pp.vystupniAdresar
  100. self.VystupFS = self.focalStatistics()
  101. self.konstanta = "-1"
  102.  
  103. def focalStatistics(self):
  104. outputFiles = []
  105. for i in range(0, len(self.VystupFILTER_LOW)):
  106. for neighbours in [1, 2, 3, 4]:
  107. for method in ["MEAN", "MAXIMUM"]:
  108. neighboursParam = "Circle " + str(neighbours) + " CELL"
  109. of = self.vystupniAdresar + "_" + str(i+1) + "_fs" + "c" + str(neighbours) + method[:2]
  110. print(self.VystupFILTER_LOW[i], of, neighboursParam, method)
  111. arcpy.gp.FocalStatistics_sa(self.VystupFILTER_LOW[i], of, neighboursParam, method, "DATA")
  112. outputFiles.append(of)
  113. return outputFiles
  114.  
  115. def times(self):
  116. outputFiles = []
  117. for i in range(0, len(self.VystupFS)):
  118. of = self.VystupFS[i] + "T"
  119. print(self.VystupFS[i], of)
  120. arcpy.Times_3d(self.VystupFS[i], self.konstanta, of)
  121. outputFiles.append(of)
  122. return outputFiles
  123.  
  124. #6 class CHMTopsDetection: vyuziti nastroju hydrologickeho modelovani,
  125. # reklasifikace a konverzi dat pro zjisteni potencialnich vrcholu stromu
  126. class CHMTopsDetection():
  127. def __init__(self):
  128. self.f_dir_drop = pp.f_dir_drop
  129. self.VystupTIMES = chmproc.times()
  130. self.FlowDir = self.flowDirection()
  131. self.FLE = self.flowLength()
  132. self.VystupniAdresarReclassify = pp.VystupniAdresarReclassify
  133. self.ReclassifiedFLE = self.reclassify()
  134. self.vystupniAdresar = pp.vystupniAdresar
  135. self.RasToPo = self.rasterToPolygon()
  136. self.Points = self.polygonToPoints()
  137. self.VystupniAdresarVYSLEDKY = pp.VystupniAdresarVYSLEDKY
  138. self.nastavPlochy = pp.nastavPlochy()
  139.  
  140. def flowDirection(self):
  141. outputFiles = []
  142. for i in range (0, len(self.VystupTIMES)):
  143. of = self.VystupTIMES[i] + "_FD"
  144. print(self.VystupTIMES[i], of, self.f_dir_drop)
  145. arcpy.gp.FlowDirection_sa(self.VystupTIMES[i], of, "FORCE", self.f_dir_drop)
  146. outputFiles.append(of)
  147. return outputFiles
  148.  
  149. def flowLength(self):
  150. outputFiles = []
  151. for i in range (0, len(self.FlowDir)):
  152. of = self.FlowDir[i][:-4] + "_FL"
  153. print(self.FlowDir[i], of)
  154. arcpy.gp.FlowLength_sa(self.FlowDir[i], of, "DOWNSTREAM", "")
  155. outputFiles.append(of)
  156. return outputFiles
  157.  
  158. def reclassify(self):
  159. outputFiles = []
  160. for i in range (0, len(self.FLE)):
  161. of = self.VystupniAdresarReclassify + "Reclass_" + self.FLE[i][-12:-3] + "_fl_l1"
  162. print(self.FLE[i], of)
  163. arcpy.Reclassify_3d(self.FLE[i], "Value", "0 1;0 999 NODATA", of, "DATA")
  164. outputFiles.append(of)
  165. return outputFiles
  166.  
  167. def rasterToPolygon(self):
  168. outputFiles = []
  169. for i in range(0, len(self.ReclassifiedFLE)):
  170. of = self.vystupniAdresar + self.ReclassifiedFLE[i][-16:-4] + "_polyg" + ".shp"
  171. print(self.ReclassifiedFLE[i], of)
  172. arcpy.RasterToPolygon_conversion(self.ReclassifiedFLE[i], of, "SIMPLIFY", "VALUE")
  173. outputFiles.append(of)
  174. return outputFiles
  175.  
  176. def polygonToPoints(self):
  177. outputFiles = []
  178. for i in range(0, len(self.RasToPo)):
  179. of = self.vystupniAdresar + self.RasToPo[i][-22:-10] + "_tops" + ".shp"
  180. print(self.RasToPo[i], of)
  181. arcpy.FeatureToPoint_management(self.RasToPo[i], of, "INSIDE")
  182. outputFiles.append(of)
  183. return outputFiles
  184.  
  185. def clipFinal(self):
  186. outputFiles = []
  187. for i in range(0, len(self.Points)):
  188. of = self.VystupniAdresarVYSLEDKY + self.Points[i][-21:-4] + "_final" + ".shp"
  189. j = i / 8
  190. print(self.Points[i], self.nastavPlochy[j], of)
  191. arcpy.Clip_analysis(self.Points[i], self.nastavPlochy[j], of, "")
  192. outputFiles.append(of)
  193. return outputFiles
  194.  
  195. #7 class FinalHeights: zjisteni vysek stromu v porostu pomoci CHM a potencialnich vrcholu stromu,
  196. # zjisteni prumeru a medianu vysek dle vsech variant nastaveni
  197. # slouceni vysledku vysek ze vsech ploch - zvlast tabulka pro kazde jednotlive nastaveni Focal Statistics
  198. # a export do excelu
  199. class FinalHeights():
  200. def __init__(self):
  201. self.TopsFinal = chmtd.clipFinal()
  202. self.VystupniAdresarVYSLEDKY = pp.VystupniAdresarVYSLEDKY
  203. self.VystupMINUS = chmpre.VystupMINUS
  204. self.HEIGHTS = self.extractValues()
  205. self.MEANS = self.meanCenter()
  206. self.MEDIANS = self.medianCenter()
  207. self.VystupniAdresarEXCEL = pp.VystupniAdresarEXCEL
  208. self.MEANS_merge = self.mergeMEANS()
  209. self.MEDIANS_merge = self.mergeMEDIANS()
  210. self.MEANS_EXCEL = self.exportMeansToExcel()
  211. self.MEDIANS_EXCEL = self.exportMediansToExcel()
  212.  
  213. def extractValues(self):
  214. outputFiles = []
  215. for i in range(0, len(self.TopsFinal)):
  216. of = self.VystupniAdresarVYSLEDKY + self.TopsFinal[i][-27:-10] + "_heights" + ".shp"
  217. j = i / 8
  218. print(self.TopsFinal[i], self.VystupMINUS[j], of)
  219. arcpy.gp.ExtractValuesToPoints_sa(self.TopsFinal[i], self.VystupMINUS[j], of, "NONE", "ALL")
  220. outputFiles.append(of)
  221. return outputFiles
  222.  
  223. def meanCenter(self):
  224. outputFiles = []
  225. for i in range(0, len(self.HEIGHTS)):
  226. of = self.VystupniAdresarVYSLEDKY + self.HEIGHTS[i][-29:-17] + "_mean_center" + ".shp"
  227. print(self.HEIGHTS[i], of)
  228. arcpy.MeanCenter_stats(self.HEIGHTS[i], of, "", "", "RASTERVALU")
  229. outputFiles.append(of)
  230. return outputFiles
  231.  
  232. def medianCenter(self):
  233. outputFiles = []
  234. for i in range(0, len(self.HEIGHTS)):
  235. of = self.VystupniAdresarVYSLEDKY + self.HEIGHTS[i][-29:-17] + "_medi_center" + ".shp"
  236. print(self.HEIGHTS[i], of)
  237. arcpy.MedianCenter_stats(self.HEIGHTS[i], of, "", "", "RASTERVALU")
  238. outputFiles.append(of)
  239. return outputFiles
  240.  
  241. def mergeMEANS(self):
  242. outputFiles = []
  243. for i in range(0, 8):
  244. of = self.VystupniAdresarVYSLEDKY + self.MEANS[i][-22:-15] + "_prumery.shp"
  245. print(self.MEANS[i::8], of)
  246. arcpy.Merge_management(self.MEANS[i::8], of, "")
  247. outputFiles.append(of)
  248. return outputFiles
  249.  
  250. def mergeMEDIANS(self):
  251. outputFiles = []
  252. for i in range(0, 8):
  253. of = self.VystupniAdresarVYSLEDKY + self.MEDIANS[i][-22:-15] + "_mediany.shp"
  254. print(self.MEDIANS[i::8], of)
  255. arcpy.Merge_management(self.MEDIANS[i::8], of, "")
  256. outputFiles.append(of)
  257. return outputFiles
  258.  
  259. def exportMeansToExcel(self):
  260. outputFiles = []
  261. for i in range(0, len(self.MEANS_merge)):
  262. of = self.VystupniAdresarEXCEL + self.MEANS_merge[i][-19:-4] + ".xls"
  263. print(self.MEANS_merge[i], of)
  264. arcpy.TableToExcel_conversion(self.MEANS_merge[i], of , "NAME", "CODE")
  265. outputFiles.append(of)
  266. return outputFiles
  267.  
  268. def exportMediansToExcel(self):
  269. outputFiles = []
  270. for i in range(0, len(self.MEDIANS_merge)):
  271. of = self.VystupniAdresarEXCEL + self.MEDIANS_merge[i][-19:-4] + ".xls"
  272. print(self.MEDIANS_merge[i], of)
  273. arcpy.TableToExcel_conversion(self.MEDIANS_merge[i], of , "NAME", "CODE")
  274. outputFiles.append(of)
  275. return outputFiles
  276.  
  277. #8 class FinalTops: zjisteni celkoveho poctu vrcholu stromu na kazde plose
  278. # slouceni vysledku (poctu vrcholu stromu) ze vsech ploch
  279. # - zvlast tabulka pro kazde jednotlive nastaveni Focal Statistics a export do excelu
  280. class FinalTops():
  281. def __init__(self):
  282. self.HEIGHTS = fh.HEIGHTS
  283. self.VystupniAdresarVYSLEDKY = pp.VystupniAdresarVYSLEDKY
  284. self.VystupTopsCount = self.TopsCount()
  285. self.mergeTopsCount = self.mergeTopsCount()
  286. self.VystupniAdresarEXCEL = pp.VystupniAdresarEXCEL
  287. self.TopsCountToExcel = self.exportTopsCountToExcel()
  288.  
  289. def TopsCount(self):
  290. outputFiles = []
  291. for i in range (0, len(self.HEIGHTS)):
  292. of = self.VystupniAdresarVYSLEDKY + self.HEIGHTS[i][-29:-19] + "_St.dbf"
  293. print(self.HEIGHTS[i], of)
  294. arcpy.Statistics_analysis(self.HEIGHTS[i], of, "FID COUNT", "")
  295. outputFiles.append(of)
  296. return outputFiles
  297.  
  298. def mergeTopsCount(self):
  299. outputFiles = []
  300. for i in range(0, 8):
  301. of = self.VystupniAdresarVYSLEDKY + self.VystupTopsCount[i][-11:-4] + "_TopsCount.dbf"
  302. print(self.VystupTopsCount[i::8], of)
  303. arcpy.Merge_management(self.VystupTopsCount[i::8], of, "")
  304. outputFiles.append(of)
  305. return outputFiles
  306.  
  307. def exportTopsCountToExcel(self):
  308. outputFiles = []
  309. for i in range(0, len(self.mergeTopsCount)):
  310. of = self.VystupniAdresarEXCEL + self.mergeTopsCount[i][-21:-4] + ".xls"
  311. print(self.mergeTopsCount[i], of)
  312. arcpy.TableToExcel_conversion(self.mergeTopsCount[i], of , "NAME", "CODE")
  313. outputFiles.append(of)
  314. return outputFiles
  315.  
  316. #9 "zavolani" vsech funkci
  317. pp = PreProcessing()
  318. pp.do_preprocessing()
  319.  
  320. chmpre = CHMPreparation()
  321. chmpre.do_CHMPreparation()
  322.  
  323. chmproc = CHMProcessing()
  324.  
  325. chmtd = CHMTopsDetection()
  326.  
  327. fh = FinalHeights()
  328.  
  329. ft = FinalTops()
  330.  
  331. print "Vsechny vypocty byly provedeny v poradku, vysledky najdete v tabulkach excel ulozenych ve slozce -uloziste_vysledky_excel- na disku C.\n"
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!