I need a script, that loads a DEM, loads a polygon with point (airports) and loads a polygon of cities. Now I have to run visibility analysis (3D analyst) from one specific airport. And then print all the cities visibile from that point in 10km perimeter and save them to new shapefile. I've coded part, but I am stuck - visibility tool always pop up error 9999: Failed to open raster dataset. Failed to execute(visibility. I am desperate and dont know what to do.
import ogrimport arcpyfrom arcpy import envfrom arcpy.da import *import gdalfrom gdalconst import *import numpyimport osimport sysfrom arcpy import envcesta_letiste = 'ArcCR500_v32/Letiste.shp'cesta_obce = 'ArcCR500_v32/ObcePolygony.shp'vystup = "kle0051_vystupy/obce_vystup.shp"vystup_mosnov = 'kle0051_vystupy/mosnov.shp'visionout = "kle0051_vystupy/visionout"env.workspace = 'ArcCR500_v32/'arcpy.env.overwriteOutput = Truedem="ArcCR500_v32/dem"if not os.path.exists(dem): print "DEM neexistuje" sys.exit()if not os.path.exists(cesta_obce): print "ObcePolygony neexistuje" sys.exit()if not os.path.exists(cesta_letiste): print "Letiste neexistuje" sys.exit()cti_obce = ogr.Open(cesta_obce)vrstva_obce = cti_obce.GetLayer(0)cti_letiste = ogr.Open(cesta_letiste)vrstva_letiste = cti_letiste.GetLayer(0)driver = ogr.GetDriverByName("ESRI Shapefile")if os.path.exists(vystup): driver.DeleteDataSource(vystup)newDataset = driver.CreateDataSource(vystup)vrstva_MSK = newDataset.CreateLayer(vystup, geom_type=ogr.wkbPolygon)vrstva_obce.SetAttributeFilter("KOD_CZNUTS='CZ080'")#zjisteni poctu obci po filtracipocet_obci = vrstva_obce.GetFeatureCount()print pocet_obcifor i in range(pocet_obci): obce_prvek= vrstva_obce.GetNextFeature() vrstva_MSK.CreateFeature(obce_prvek) exampleFeature = vrstva_MSK.GetFeature(0)newDataset.Destroy()driver = ogr.GetDriverByName("ESRI Shapefile")if os.path.exists(vystup_mosnov): driver.DeleteDataSource(vystup_mosnov)dataset2=driver.CreateDataSource(vystup_mosnov)mosnov_bod=dataset2.CreateLayer(vystup_mosnov, geom_type = ogr.wkbPoint)vrstva_letiste.SetAttributeFilter("NAZEV_ASCI='Ostrava/Mosnov'")mosnov = ogr.Geometry(ogr.wkbPoint)mosnov_prvek = vrstva_letiste.GetNextFeature()mosnov_letiste = mosnov_prvek.GetGeometryRef()mosnov_bod.CreateFeature(mosnov_prvek)arcpy.CheckOutExtension("3D")#vytvoreni bufferu pro mosnov (okoli 10km)mosnov_buffer = mosnov_letiste.Buffer(10000)algOutput = "kle0051_vystupy/algOutput1"analysisType = "FREQUENCY"nonVisibleValue = "ZERO"zFactor = 1useEarthCurvature = "CURVED_EARTH"refractivityCoefficient = 0.13surfaceOffset = 500observerElevation = 2000observerOffset = 500innerRadius = 20000outerRadius = 100000horizStartAngle = 45horizEndAngle = 215vertUpperAngle = 5vertLowerAngle = -5arcpy.Visibility_3d(dem, vystup_mosnov, visionout, algOutput, analysisType, nonVisibleValue, zFactor, useEarthCurvature, refractivityCoefficient, surfaceOffset, observerElevation, observerOffset, innerRadius, outerRadius, horizStartAngle, horizEndAngle, vertUpperAngle, vertLowerAngle)
أكثر...
import ogrimport arcpyfrom arcpy import envfrom arcpy.da import *import gdalfrom gdalconst import *import numpyimport osimport sysfrom arcpy import envcesta_letiste = 'ArcCR500_v32/Letiste.shp'cesta_obce = 'ArcCR500_v32/ObcePolygony.shp'vystup = "kle0051_vystupy/obce_vystup.shp"vystup_mosnov = 'kle0051_vystupy/mosnov.shp'visionout = "kle0051_vystupy/visionout"env.workspace = 'ArcCR500_v32/'arcpy.env.overwriteOutput = Truedem="ArcCR500_v32/dem"if not os.path.exists(dem): print "DEM neexistuje" sys.exit()if not os.path.exists(cesta_obce): print "ObcePolygony neexistuje" sys.exit()if not os.path.exists(cesta_letiste): print "Letiste neexistuje" sys.exit()cti_obce = ogr.Open(cesta_obce)vrstva_obce = cti_obce.GetLayer(0)cti_letiste = ogr.Open(cesta_letiste)vrstva_letiste = cti_letiste.GetLayer(0)driver = ogr.GetDriverByName("ESRI Shapefile")if os.path.exists(vystup): driver.DeleteDataSource(vystup)newDataset = driver.CreateDataSource(vystup)vrstva_MSK = newDataset.CreateLayer(vystup, geom_type=ogr.wkbPolygon)vrstva_obce.SetAttributeFilter("KOD_CZNUTS='CZ080'")#zjisteni poctu obci po filtracipocet_obci = vrstva_obce.GetFeatureCount()print pocet_obcifor i in range(pocet_obci): obce_prvek= vrstva_obce.GetNextFeature() vrstva_MSK.CreateFeature(obce_prvek) exampleFeature = vrstva_MSK.GetFeature(0)newDataset.Destroy()driver = ogr.GetDriverByName("ESRI Shapefile")if os.path.exists(vystup_mosnov): driver.DeleteDataSource(vystup_mosnov)dataset2=driver.CreateDataSource(vystup_mosnov)mosnov_bod=dataset2.CreateLayer(vystup_mosnov, geom_type = ogr.wkbPoint)vrstva_letiste.SetAttributeFilter("NAZEV_ASCI='Ostrava/Mosnov'")mosnov = ogr.Geometry(ogr.wkbPoint)mosnov_prvek = vrstva_letiste.GetNextFeature()mosnov_letiste = mosnov_prvek.GetGeometryRef()mosnov_bod.CreateFeature(mosnov_prvek)arcpy.CheckOutExtension("3D")#vytvoreni bufferu pro mosnov (okoli 10km)mosnov_buffer = mosnov_letiste.Buffer(10000)algOutput = "kle0051_vystupy/algOutput1"analysisType = "FREQUENCY"nonVisibleValue = "ZERO"zFactor = 1useEarthCurvature = "CURVED_EARTH"refractivityCoefficient = 0.13surfaceOffset = 500observerElevation = 2000observerOffset = 500innerRadius = 20000outerRadius = 100000horizStartAngle = 45horizEndAngle = 215vertUpperAngle = 5vertLowerAngle = -5arcpy.Visibility_3d(dem, vystup_mosnov, visionout, algOutput, analysisType, nonVisibleValue, zFactor, useEarthCurvature, refractivityCoefficient, surfaceOffset, observerElevation, observerOffset, innerRadius, outerRadius, horizStartAngle, horizEndAngle, vertUpperAngle, vertLowerAngle)
أكثر...