ArcPy to Rasterio: A Dated Technical Migration Note

import arcpy
from arcpy.sa import *

# Set up environment
arcpy.env.workspace = "C:/data"
arcpy.CheckOutExtension("Spatial")

# Read raster - returns Raster object
dem = Raster("elevation.tif")

# Access properties
print(f"Cell size: {dem.meanCellWidth}")
print(f"Extent: {dem.extent}")
print(f"Spatial ref: {dem.spatialReference.name}")

import rasterio
import numpy as np

# Open raster - context manager handles cleanup
with rasterio.open("elevation.tif") as src:
    # Read as NumPy array (band 1)
    dem = src.read(1)

    # Access properties via profile
    print(f"Cell size: {src.res}")
    print(f"Bounds: {src.bounds}")
    print(f"CRS: {src.crs}")

    # dem is now a standard NumPy array
    print(f"Shape: {dem.shape}")
    print(f"Dtype: {dem.dtype}")

from arcpy.sa import Slope

# Calculate slope in degrees
dem = Raster("elevation.tif")
slope = Slope(dem, "DEGREE")
slope.save("slope.tif")

import rasterio
import numpy as np

with rasterio.open("elevation.tif") as src:
    dem = src.read(1).astype(float)
    cellsize = src.res[0]  # Assumes square cells

    # Calculate gradient using NumPy
    dy, dx = np.gradient(dem, cellsize)
    slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
    slope_deg = np.degrees(slope_rad)

    # Write output
    profile = src.profile.copy()
    profile.update(dtype=rasterio.float32)
    with rasterio.open("slope.tif", "w", **profile) as dst:
        dst.write(slope_deg.astype(np.float32), 1)

from arcpy.sa import Con, Reclassify, RemapRange

# Conditional: Set values < 0 to NoData
dem = Raster("elevation.tif")
dem_positive = Con(dem >= 0, dem)

# Reclassify elevation to categories
reclass = Reclassify(dem, "VALUE",
    RemapRange([[0, 100, 1],
                [100, 500, 2],
                [500, 1000, 3],
                [1000, 9999, 4]])
)
reclass.save("elevation_classes.tif")

import rasterio
import numpy as np

with rasterio.open("elevation.tif") as src:
    dem = src.read(1).astype(float)
    nodata = src.nodata

    # Conditional: Set negative values to NoData
    dem[dem < 0] = np.nan

    # Reclassify using np.digitize or np.select
    conditions = [
        (dem >= 0) & (dem < 100),
        (dem >= 100) & (dem < 500),
        (dem >= 500) & (dem < 1000),
        dem >= 1000
    ]
    classes = np.select(conditions, [1, 2, 3, 4], default=0)

    # Write output
    profile = src.profile.copy()
    profile.update(dtype=rasterio.uint8, nodata=0)
    with rasterio.open("elevation_classes.tif", "w", **profile) as dst:
        dst.write(classes.astype(np.uint8), 1)

import rasterio

# Standard GeoTIFF write
with rasterio.open(
    "output.tif", "w",
    driver="GTiff",
    height=data.shape[0],
    width=data.shape[1],
    count=1,
    dtype=data.dtype,
    crs=src.crs,
    transform=src.transform
) as dst:
    dst.write(data, 1)

# Problem: No internal tiling
# Problem: No overviews
# Problem: Must download entire file to read

import rasterio
from rasterio.enums import Resampling

# Cloud-Optimized GeoTIFF with COG driver
with rasterio.open(
    "output_cog.tif", "w",
    driver="COG",
    height=data.shape[0],
    width=data.shape[1],
    count=1,
    dtype=data.dtype,
    crs=src.crs,
    transform=src.transform,
    compress="deflate",        # Compression
    blocksize=512,             # Internal tile size
    overview_resampling="average"
) as dst:
    dst.write(data, 1)

# Or use rio-cogeo CLI:
# rio cogeo create input.tif output_cog.tif

# DON'T: Hardcode metadata
with rasterio.open(
    "output.tif", "w",
    driver="GTiff",
    height=1000,
    width=1000,
    count=1,
    dtype="float32",
    crs="EPSG:4326",  # Might not match input!
    transform=Affine(0.1, 0, 0, 0, -0.1, 0)  # Wrong!
) as dst:
    dst.write(data, 1)

# Problems:
# - CRS might not match input
# - Transform is completely wrong
# - No compression (huge files)
# - nodata not set

# DO: Copy and update profile
with rasterio.open("input.tif") as src:
    data = src.read(1)
    profile = src.profile.copy()

    # Process data
    result = some_processing(data)

    # Update only what changed
    profile.update(
        driver="COG",           # Cloud-optimized
        dtype=result.dtype,     # Match output type
        compress="deflate",     # Compression
        nodata=-9999           # Explicit nodata
    )

    with rasterio.open("output.tif", "w", **profile) as dst:
        dst.write(result, 1)