Python rozhraní Esri ArcGIS Pro - arcpy¶
Ukázka využití arcpy vychází z předmětu GIS 2.
Vstupní data byla získána z EU katalogu: NUTS - Nomenclature of territorial units for statistics classification (vydavatel: Eurostat).
Vektorová data¶
Základní informace o datových zdrojích poskytuje Describe. V případě níže jde o Feature Class a Table.
nuts_path = r"S:\K155\Public\155UZPR\01\NUTS_RG_20M_2021_3035.shp"
nuts = arcpy.Describe(nuts_path)
nuts.shapeType
'Polygon'
Atributy a jejich vlastnosti jsou dostupné pomocí fields
(viz Table).
len(nuts.fields)
11
for field in nuts.fields:
print(field.name, field.type)
FID OID Shape Geometry NUTS_ID String LEVL_CODE SmallInteger CNTR_CODE String NAME_LATN String NUTS_NAME String MOUNT_TYPE SmallInteger URBN_TYPE SmallInteger COAST_TYPE SmallInteger FID_1 String
Záznamy v atributové tabulce jsou dostupné pomocí tzv. SearchCursor
.
rows = arcpy.SearchCursor(nuts_path)
for row in rows:
print(row.CNTR_CODE)
break
FR
count = {}
rows = arcpy.SearchCursor(nuts_path)
for row in rows:
if row.LEVL_CODE not in count:
count[row.LEVL_CODE] = 0
count[row.LEVL_CODE] += 1
print(count)
{0: 37, 1: 125, 2: 334, 3: 1514}
Získejme přehled všech kódů jednotlivých států.
codes = []
rows = arcpy.SearchCursor(nuts_path)
count = 0
for row in rows:
if row.CNTR_CODE not in codes:
codes.append(row.CNTR_CODE)
print(sorted(codes))
['AL', 'AT', 'BE', 'BG', 'CH', 'CY', 'CZ', 'DE', 'DK', 'EE', 'EL', 'ES', 'FI', 'FR', 'HR', 'HU', 'IE', 'IS', 'IT', 'LI', 'LT', 'LU', 'LV', 'ME', 'MK', 'MT', 'NL', 'NO', 'PL', 'PT', 'RO', 'RS', 'SE', 'SI', 'SK', 'TR', 'UK']
Geometrie je přístupná v podobě atributu SHAPE
.
code = 'CZ'
area = 0
rows = arcpy.SearchCursor(nuts_path, f"CNTR_CODE = '{code}' AND LEVL_CODE=0")
for row in rows:
area += row.SHAPE.area
print(f"{code} area in km2: {area/1e6:.2f}")
CZ area in km2: 78896.87
ArcPy umožňuje spoustět jakýkoliv nástroj, který je možný spustit v rámci grafického uživatelského rozhraní ArcGIS Pro. V příkladu níže Select. V nápovědě k nástroji najdete příklad použití v Pythonu.
arcpy.env.overwriteOutput = True
for code in codes[:3]:
print("{}...".format(code))
arcpy.analysis.Select(nuts_path, f"nuts_{code}", f"CNTR_CODE = '{code}' AND LEVL_CODE = 2")
FR... HR... HU...
for code in codes[:3]:
print("{}...".format(code))
arcpy.management.Dissolve(f"nuts_{code}", f"nuts0_{code}", "CNTR_CODE")
FR... HR... HU...
Rastrová data¶
Základní informace o rastrových datech opět poskytuje Describe. V tomto případě jde o V případě níže jde o Raster Dataset.
dem_path = r"S:\K155\Public\155UZPR\01\eu_dem_v11_E40N30.tif"
dem = arcpy.Describe(dem_path)
dem.width, dem.height
(40000, 40000)
Zájmové území vybereme na základě NUTS jednotky a vrstvu digitálního modelu terénu na základě toho ořežeme pomocí nástroje Clip.
nuts_name = 'Ústecký kraj'
aoi_path = "nuts_aoi"
arcpy.analysis.Select(nuts_path, aoi_path, f"NUTS_NAME = '{nuts_name}'")
dem_aoi_path = "dmt_aoi"
arcpy.management.Clip(dem_path, out_raster=dem_aoi_path, in_template_dataset=aoi_path,
clipping_geometry="ClippingGeometry")
dem_aoi = arcpy.Describe(dem_aoi_path)
dem_aoi.width, dem_aoi.height
(4618, 4426)
Na základě rastrových dat můžeme vytvořit numpy pole a dále pro výpočet používat knihovnu Numpy.
dem_aoi_data = arcpy.Raster(dem_aoi_path)
array = arcpy.RasterToNumPyArray(dem_aoi_data)
array.shape, array.dtype
((4426, 4618), dtype('float32'))
Získejme minimální a maximální nadmořskou výšku pro dané zájmové území.
array.min(), array.max()
(-3.402823e+38, 1228.7755)
import numpy
array = arcpy.RasterToNumPyArray(dem_aoi_data, nodata_to_value=numpy.nan)
numpy.nanmin(array), numpy.nanmax(array)
(49.590366, 1228.7755)
Pomocí numpy můžeme provést např. reklasifikaci hodnot.
array[array<=200] = 0
array[array>200] = 1
numpy.unique(array, return_counts=True) # numpy >= 1.24 equal_nan=True
(array([ 0., 1., nan, ..., nan, nan, nan], dtype=float32), array([ 913759, 7596521, 1, ..., 1, 1, 1], dtype=int64))
Výsledek reklasifikace rastrových hodnot zapíšeme na disk.
arcpy.env.outputCoordinateSystem = dem_aoi_path
arcpy.env.cellSize = dem_aoi_path
no_data_value = 9999
array[numpy.isnan(array)] = no_data_value
output = arcpy.NumPyArrayToRaster(array.astype('int16'), arcpy.Point(dem_aoi.extent.XMin, dem_aoi.extent.YMin), dem_aoi.meanCellWidth, dem_aoi.meanCellHeight, value_to_nodata=no_data_value)
output.save("dmt200")
Výsledek výpočtu lze zobrazit přímo v notebooku