In this tutorial, we will learn about Arcpy, which is a popular Python library used in geographic information system. Here we will explore geoprocessing tools by importing Arcpy also we will learn some indices such as NDVI, NDBI, NDWI, NDMI.
N.B. Currently ArcPy is not readily supported in Google Colab. So please use these codes in the ArcGIS Python Environment.
Creating a Folder and Add a shapefile¶
1. Creating a Folder¶
# Specify the path of the parent directory
parent_path = "e:/"
# Specify the name of the new folder
Folder_name = "PythonExercise"
# Create the new folder
arcpy.CreateFolder_management(parent_path,Folder_name)
2. Add Shapefile as a Layer¶
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the path to the shapefile
shapefile = "E:\Files_E\GIS Data Documentation\BD Data"
# Add shapefile as a layer
arcpy.MakeFeatureLayer_management(shapefile,"Sylhet")
Geoprocessing Tool¶
1. Buffer Analysis¶
A buffer refers to a zone or area that is created around a specific geographic feature, such as a point, line, or polygon. This buffer zone is defined by a specified distance or radius and is used for various analytical and visualization purposes.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the path to the input feature
input_feature = "E:\Files_E\GIS Data Documentation\BD Data\Sylhet.shp"
# Specify the output location and name for the buffer
output_feature = "E:\Files_E\GIS Data Documentation\BD Data\Sylhet_Buffer.shp"
# Specify the buffer distance
buffer_distance = "10 Kilometers"
# Perform the buffer analysis
arcpy.Buffer_analysis(input_feature,output_feature,buffer_distance)
2. Clip¶
Clip refers to a spatial operation that involves cutting or extracting a portion of a geographic dataset or shapefile based on the boundaries of another geographic dataset or shapefile. The process essentially clips or trims the original dataset to retain only the areas that fall within the boundary of the second dataset.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the path to the input feature class
input_feature = "E:\Files_E\GIS Data Documentation\BD Data\BD_districts.shp"
# Define the path to the clip feature
clip_feature = "E:\Files_E\GIS Data Documentation\BD Data\Sylhet.shp"
# Specify the output location and name for the clipped feature class
output_feature = "E:\Files_E\GIS Data Documentation\BD Data\Sylhet_Clip.shp"
# Perform the clip analysis
arcpy.Clip_analysis(input_feature,clip_feature,output_feature)
3. Intersect¶
Intersect is a spatial operation that involves finding and retaining only the areas or features that are common to two or more geographic datasets. The intersect operation creates a new dataset that includes only the spatially overlapping portions of the input datasets.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the paths to the input feature classes
input_feature_classes = ["E:\Files_E\GIS Data Documentation\BD Data\BD_districts.shp", "E:\Files_E\GIS Data Documentation\BD Data\Sylhet.shp"]
# Specify the output location and name for the intersected feature class
output_feature_class = "E:\Files_E\GIS Data Documentation\BD Data\Sylhet_Intersect.shp"
# Perform the intersect analysis
arcpy.Intersect_analysis(input_feature_classes, output_feature_class)
4. Union¶
Union is a spatial operation that combines the geometry of two or more geographic datasets to create a new dataset that includes all the features from the input layers.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the paths to the input feature classes
input_feature_classes = ["E:\Files_E\GIS Data Documentation\BD Data\Sylhet.shp", "E:\Files_E\GIS Data Documentation\BD Data\Sunamganj.shp"]
# Specify the output location and name for the unioned feature class
output_feature_class = "E:\Files_E\GIS Data Documentation\BD Data\Sunamganj_Sylhet_Union.shp"
# Perform the union analysis
arcpy.Union_analysis(input_feature_classes, output_feature_class)
5. Dissolve¶
Dissolve is a spatial operation that merges or aggregates features within a single layer based on a common attribute or a shared characteristic. The dissolve operation simplifies the geometry of a dataset by combining adjacent or overlapping features that have the same attribute value, resulting in a new dataset with fewer boundaries and a reduced number of features.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the path to the input feature class
input_feature_class = "E:\Files_E\GIS Data Documentation\BD Data\BD_districts.shp"
# Specify the output location and name for the dissolved feature class
output_feature_class = "E:\test\Sylhet_Division_dissolve.shp"
# Specify the field(s) to dissolve
dissolve_fields = ["DistrictName"]
# Perform the dissolve analysis
arcpy.Dissolve_management(input_feature_class, output_feature_class, dissolve_fields)
6. Merge¶
Merge is a data processing operation that combines two or more separate geographic datasets into a single dataset. The merge operation brings together the attributes and geometries from the input layers to create a unified output layer with all the features and information from the original datasets.
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"E:\Files_E\GIS Data Documentation"
arcpy.env.overwriteOutput = True
# Define the paths to the input feature classes
input_feature_classes = ["E:\Files_E\GIS Data Documentation\BD Data\Sylhet.shp", "E:\Files_E\GIS Data Documentation\BD Data\Sunamganj.shp"]
# Specify the output location and name for the merged feature class
output_feature_class = "E:\Files_E\GIS Data Documentation\BD Data\Sunamganj_Sylhet_Merge.shp"
# Perform the merge analysis
arcpy.Merge_management(input_feature_classes, output_feature_class)
7. Raster to Polygon¶
import arcpy
# Set up the workspace and overwrite output
arcpy.env.workspace = r"D:\Files_D\Python Series\Conversion"
arcpy.env.overwriteOutput = True
# Define the path to the input raster dataset
input_raster = r"D:\Files_D\Watershed Delination\basin\basin.tif"
# Define the path and name for the output polygon feature class
output_polygon = r"D:\Files_D\Python Series\Conversion\Polygon\Raster_to_Polygon"
# Perform the Raster to Polygon conversion
arcpy.RasterToPolygon_conversion(input_raster, output_polygon, simplify="NO_SIMPLIFY")
Geodatabase¶
A geodatabase is a container or data storage format used to store, manage, and organize geographic data.
1. Creating a Geodatabase¶
import arcpy
# Set up the workspace
workspace = r"E:\Files_E\Arcpy\Geodatabase"
arcpy.env.workspace = workspace
# Define the path and name for the new geodatabase
gdb_name = "Geodatabase.gdb"
gdb_path = arcpy.CreateFileGDB_management(workspace, gdb_name)
print("Geodatabase created at: {}".format(gdb_path))
2. Add shapefile to the Geodatabase¶
import arcpy
# Set up the workspace
workspace = r"E:\Files_E\Arcpy\Geodatabase"
arcpy.env.workspace = workspace
# Define the path to the input shapefile
input_shapefile = r"E:\Files_E\GIS Data Documentation\BD Data\BD_admin_boundary.shp"
# Define the name for the imported feature class in the geodatabase
output_feature_class = "Feature"
# Add the shapefile to the geodatabase
arcpy.CopyFeatures_management(input_shapefile, output_feature_class)
NDWI¶
Normalized Difference Water Index (NDWI) is used to assess the presence and extent of water bodies in satellite imagery. The index is especially useful for monitoring changes in water bodies, such as lakes, rivers, and reservoirs.
The formula to calculate NDWI is as follows:
NDWI = (Green - NIR) / (Green + NIR)
Where:
Green is the reflectance in the green band.
NIR stands for Near-Infrared reflectance.
The values of NDWI typically range from -1 to +1, where negative values indicate non-water areas, and positive values indicate water presence. The higher the positive value, the more significant the presence of water in the observed area.
# Import the necessary Arcpy modules
import arcpy
from arcpy.sa import *
# Set up the workspace and overwrite output
arcpy.env.workspace = r"D:\Files_D\Python Series\NDWI"
arcpy.env.overwriteOutput = True
# Define the paths to input bands (NIR and green) for the NDWI calculation
nir_band = r"D:\Files_D\Python Series\NDWI\Raw Data\1999\LT05_L1TP_137043_19991116_20161217_01_T1_B4.TIF"
green_band = r"D:\Files_D\Python Series\NDWI\Raw Data\1999\LT05_L1TP_137043_19991116_20161217_01_T1_B2.TIF"
# Read the raster bands as raster objects
nir_raster = arcpy.Raster(nir_band)
green_raster = arcpy.Raster(green_band)
# Convert to floating-point to avoid integer division issues
nir_float = Float(nir_raster)
green_float = Float(green_raster)
# Calculate the NDWI
ndwi = (green_float - nir_float) / (green_float + nir_float)
# Save the NDWI result to a new raster file
output_path = r"D:\Files_D\Python Series\NDWI\Output\NDWI.tif"
ndwi.save(output_path)
NDBI¶
Normalized Difference Built-up Index (NDBI) is used in satellite image analysis to identify and map built-up or urban areas. NDBI helps to distinguish between built-up areas and natural features like vegetation and bare soil based on their different spectral responses.
The formula to calculate NDBI is as follows:
NDBI = (SWIR - NIR) / (SWIR + NIR)
Where:
SWIR stands for Shortwave Infrared reflectance.
NIR stands for Near-Infrared reflectance.
The index takes advantage of the fact that built-up or urban areas typically have a high reflectance in the SWIR region due to their construction materials (e.g., concrete, asphalt), while vegetation has a higher reflectance in the NIR region. As a result, NDBI can highlight urban areas that have a strong SWIR response and a weak NIR response.
The values of NDBI usually range from -1 to +1, where negative values indicate natural features (e.g., vegetation) and positive values indicate built-up areas. The higher the positive value, the more significant the presence of built-up structures.
# Import the necessary Arcpy modules
import arcpy
from arcpy.sa import *
# Set up the workspace and overwrite output
arcpy.env.workspace = r"D:\Files_D\Python Series\NDBI"
arcpy.env.overwriteOutput = True
# Define the paths to input bands (near-infrared and shortwave infrared) for the NDBI calculation
nir_band = r"D:\Files_D\Python Series\NDBI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B5.TIF"
swir_band = r"D:\Files_D\Python Series\NDBI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B6.TIF"
# Read the raster bands as raster objects
nir_raster = arcpy.Raster(nir_band)
swir_raster = arcpy.Raster(swir_band)
# Convert to floating-point to avoid integer division issues
nir_float = Float(nir_raster)
swir_float = Float(swir_raster)
# Calculate the NDBI
ndbi = (swir_float - nir_float) / (swir_float + nir_float)
# Save the NDBI result to a new raster file
output_path = r"D:\Files_D\Python Series\NDBI\Output\NDBI.tif"
ndbi.save(output_path)
NDMI¶
Normalized Difference Moisture Index (NDMI) is used to determine vegetation water content. NDMI quantifies the difference between the near-infrared (NIR) and shortwave infrared (SWIR) reflectance of vegetation.
The formula to calculate NDMI is as follows:
NDMI = (NIR - SWIR) / (NIR + SWIR)
Where:
NIR stands for Near-Infrared reflectance.
SWIR stands for Shortwave Infrared reflectance.
The NDMI values typically range from -1 to +1, where negative values represent dry or sparse vegetation, and positive values indicate healthier and more moisture-rich vegetation. Vegetation with higher moisture content will have higher positive NDMI values.
# Import the necessary Arcpy modules
import arcpy
from arcpy.sa import *
# Set up the workspace and overwrite output
arcpy.env.workspace = r"D:\Files_D\Python Series\NDMI"
arcpy.env.overwriteOutput = True
# Define the paths to input bands (shortwave infrared 1 and near-infrared) for the NDMI calculation
swir1_band = r"D:\Files_D\Python Series\NDMI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B6.TIF"
nir_band = r"D:\Files_D\Python Series\NDMI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B5.TIF"
# Read the raster bands as raster objects
swir1_raster = arcpy.Raster(swir1_band)
nir_raster = arcpy.Raster(nir_band)
# Convert to floating-point to avoid integer division issues
swir1_float = Float(swir1_raster)
nir_float = Float(nir_raster)
# Calculate the NDMI
ndmi = (nir_float - swir1_float) / (nir_float + swir1_float)
# Save the NDMI result to a new raster file
output_path = r"D:\Files_D\Python Series\NDMI\Output\NDMI.tif"
ndmi.save(output_path)
NDVI¶
Normalized Difference Vegetation Index (NDVI) is to assess the health, density, and vigor of vegetation cover on the Earth's surface. NDVI quantifies the difference between the near-infrared (NIR) and red reflectance of vegetation.
The formula to calculate NDVI is as follows:
NDVI = (NIR - Red) / (NIR + Red)
Where:
NIR stands for Near-Infrared reflectance.
Red stands for Red reflectance.
The NDVI values typically range from -1 to +1, where negative values represent non-vegetated areas (like water bodies or barren land), values close to zero represent sparse or stressed vegetation, and positive values indicate healthy and dense vegetation. Lush and healthy vegetation will have higher positive NDVI values.
# Import the necessary Arcpy modules
import arcpy
from arcpy.sa import *
# Set up the workspace and overwrite output
arcpy.env.workspace = r"D:\Files_D\Python Series\NDVI"
arcpy.env.overwriteOutput = True
# Define the paths to input bands (near-infrared and red bands) for the NDVI calculation
nir_band = r"D:\Files_D\Python Series\NDVI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B5.TIF"
red_band = r"D:\Files_D\Python Series\NDVI\Raw Data\2023\LC09_L2SP_136043_20230119_20230313_02_T1_SR_B4.TIF"
# Read the raster bands as raster objects
nir_raster = arcpy.Raster(nir_band)
red_raster = arcpy.Raster(red_band)
# Convert to floating-point to avoid integer division issues
nir_float = Float(nir_raster)
red_float = Float(red_raster)
# Calculate the NDVI
ndvi = (nir_float - red_float) / (nir_float + red_float)
# Save the NDVI result to a new raster file
output_path = r"D:\Files_D\Python Series\NDVI\Output\NDVI.tif"
ndvi.save(output_path)
