4.2.2 Raster data processing¶

Clip DEM by selected NUTS region.

In [31]:

arcpy.env.workspace = r"C:\Users\martin\Documents\isdp_04" # change path here
nuts_name = "Potenza"
aoi = f"nuts_{nuts_name}.shp"
arcpy.analysis.Select(nuts_path, aoi, f"NUTS_NAME = '{nuts_name}' AND LEVL_CODE = 3")

dmt =  r"S:\K155\Public\155ISDP\01\eu_dem_v11_E40N10.TIF"
dmt_aoi = "dem_aoi.TIF"
arcpy.management.Clip(dmt, out_raster=dmt_aoi, in_template_dataset=aoi,
                      clipping_geometry="ClippingGeometry")

arcpy.env.workspace = r"C:\Users\martin\Documents\isdp_04" # change path here nuts_name = "Potenza" aoi = f"nuts_{nuts_name}.shp" arcpy.analysis.Select(nuts_path, aoi, f"NUTS_NAME = '{nuts_name}' AND LEVL_CODE = 3") dmt = r"S:\K155\Public\155ISDP\01\eu_dem_v11_E40N10.TIF" dmt_aoi = "dem_aoi.TIF" arcpy.management.Clip(dmt, out_raster=dmt_aoi, in_template_dataset=aoi, clipping_geometry="ClippingGeometry")

Out[31]:

Messages

Start Time: Monday, March 10, 2025 11:49:42 AM
Building Pyramids...
Succeeded at Monday, March 10, 2025 11:51:27 AM (Elapsed Time: 1 minutes 45 seconds)

Basic information about data sources provides Describe (Dataset properties).

In [32]:

data = arcpy.Describe(dmt_aoi)
data.bandCount

data = arcpy.Describe(dmt_aoi) data.bandCount

Out[32]:

1

Load raster data and check the metadata

In [33]:

# create a raster layer
raster = arcpy.Raster(dmt_aoi)

raster.width, raster.height

# create a raster layer raster = arcpy.Raster(dmt_aoi) raster.width, raster.height

Out[33]:

(3861, 4451)

Create a numpy array from the raster

In [34]:

array = arcpy.RasterToNumPyArray(raster)
type(array)

array = arcpy.RasterToNumPyArray(raster) type(array)

Out[34]:

numpy.ndarray

Check if it went well

In [35]:

array.shape, array.dtype

array.shape, array.dtype

Out[35]:

((4451, 3861), dtype('float32'))

Check minimum and maximum value

In [36]:

array.min(), array.max()

array.min(), array.max()

Out[36]:

(-3.402823e+38, 2162.334)

In [37]:

array = arcpy.RasterToNumPyArray(raster, nodata_to_value=numpy.nan)
numpy.nanmin(array), numpy.nanmax(array)

array = arcpy.RasterToNumPyArray(raster, nodata_to_value=numpy.nan) numpy.nanmin(array), numpy.nanmax(array)

Out[37]:

(-0.78258896, 2162.334)

Perform the map algebra

In [38]:

array[array<=200] = 0
array[array>200] = 1
numpy.unique(array, return_counts=True)

array[array<=200] = 0 array[array>200] = 1 numpy.unique(array, return_counts=True)

Out[38]:

(array([ 0.,  1., nan], dtype=float32),
 array([  43117, 9314997, 7827197], dtype=int64))

Write the array into a raster

In [40]:

arcpy.env.outputCoordinateSystem = dmt_aoi
arcpy.env.cellSize = dmt_aoi

no_data_value = 9999
array[numpy.isnan(array)] = no_data_value
output = arcpy.NumPyArrayToRaster(array.astype('int16'), arcpy.Point(raster.extent.XMin, raster.extent.YMin),  raster.meanCellWidth, raster.meanCellHeight, value_to_nodata=no_data_value)
output.save("dmt_200.tif")

arcpy.env.outputCoordinateSystem = dmt_aoi arcpy.env.cellSize = dmt_aoi no_data_value = 9999 array[numpy.isnan(array)] = no_data_value output = arcpy.NumPyArrayToRaster(array.astype('int16'), arcpy.Point(raster.extent.XMin, raster.extent.YMin), raster.meanCellWidth, raster.meanCellHeight, value_to_nodata=no_data_value) output.save("dmt_200.tif")

In [41]:

output

output

Out[41]: