Calculate Impervious Surfaces from Spectral Imagery

Introduction and objective

Ground surfaces that are impenetrable to water can cause serious environmental problems, including flooding and contaminated runoff. Many local government agencies use impervious surface calculations to compute the storm water bill for properties. That depends on how much impervious area is there in each property.

Map highlighting impervious surface areas: driveways, private walks, roof tops, and parking lots [1].

In this notebook, we’ll use a high resolution land cover map obtained from Chesapeake Conservancy to determine which parts of the ground are pervious and impervious.

Impervious surfaces are generally human-made: buildings, roads, parking lots, brick, or asphalt. Pervious surfaces include vegetation, water bodies, and bare soil.

We will reclassify those land-use types into either impervious or pervious surfaces to derive the stormwater bill for each property. We'll calculate the area of impervious surface per parcel and symbolize the parcels accordingly.

Necessary Imports

Input
import pandas as pd
from datetime import datetime as dt

import arcgis
import arcpy
from arcgis.gis import GIS
from arcgis.raster.functions import RFT           

gis = GIS('https://pythonapi.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')

Get data for analysis

For this analysis we need the following datasets:

  • Land Cover of Kent County, Delaware
  • Kent County parcel feature layer
  • Raster Function Template to extract necessary class labels from the land cover layer

In the cells below, we access these datasets from the GIS.

Input
rft_to_reclassify = gis.content.get('eb9fe459cc26486eb4a733b2a2def44f')
rft_to_reclassify
Output
RFT_Reclassify_RasterFunction
RFT_Reclassify_RasterFunctionRaster function template by api_data_owner
Last Modified: June 07, 2021
0 comments, 104 views
Input
classified_raster = gis.content.get('7e00802c927a40cbbdc36cefdb2928e4')
classified_raster
Output
Land Cover of Kent County Delaware
KENT_10001_us8bit_tifImagery Layer by api_data_owner
Last Modified: June 08, 2021
0 comments, 0 views
Input
classified_raster.layers[0]
Output

Twelve land cover classes were mapped:

  1. Water
  2. Emergent Wetlands
  3. Tree Canopy
  4. Scrub/Shrub
  5. Low Vegetation
  6. Barren
  7. Structures
  8. Other Impervious Surfaces
  9. Roads
  10. Tree Canopy over Structures
  11. Tree Canopy over Other Impervious Surfaces
  12. Tree Canopy over Roads

The complete class definitions and standards can be viewed at the link below.http://goo.gl/THacgg.

Input
parcels = gis.content.get('1a7e4a24372048048a668c8a176932bf')
Input
parcels
Output
kent_county_parcels
Feature Layer Collection by api_data_owner
Last Modified: February 28, 2020
0 comments, 20 views
Input
parcel_lyr = parcels.layers[0]

Classify pixel values into pervious and impervious class

Create the RFT object from raster function template portal item. This object serves as a python function corresponding to the Raster Function Template.

Input
rft_to_reclassify_ob = RFT(rft_to_reclassify)

Once an RFT class object has been created, a single question mark before, or after the RFT object will show help relative to it. The help document would display the parameters that were marked as public by the author of the RFT.

Input
rft_to_reclassify_ob?
"""
RFT to classify pixel values to pervious and impervious class

Parameters
----------

Type : 3
Rows : 5

Returns
-------
Imagery Layer, on which the function chain is applied 
"""

We can visualize the function chain and the help document specific to this RFT inorder to understand the input parameters and the functions involved.

To view the raster function chain visually, we install graphviz Python library.

Input
# ! conda install graphviz -y
Input
# rft_to_reclassify_ob.draw_graph()
Output

Provide the values for the input parameters in the RFT object to create the Imagery Layer.

Input
output = rft_to_reclassify_ob(Raster=classified_raster.layers[0])
Input
persisted_output = output.save('generate_reclassified_raster' + str(dt.now().microsecond))
Input
impervious_lyr = persisted_output.layers[0]
Input
impervious_lyr
Output

Tabulate area of impervious surfaces

We'll determine the area of impervious surfaces within each parcel of land in the neighborhood. We'll first calculate the area and store the results in a stand-alone table. Then, we'll join the table to the Parcels layer.

Input
out_path = '/arcgis/home'
Input
arcpy.CreateFileGDB_management(out_path, "output_table.gdb")

Impervious field shows the area (in feet) of impervious surfaces per parcel, while Pervious shows the area of pervious surfaces. We will calculate the impervious area with the name field which denotes Parcel_ID.

Input
tabulate_op = arcpy.sa.TabulateArea(in_zone_data=parcel_lyr.url, 
                      zone_field="name", 
                      in_class_data=impervious_lyr.url, 
                      class_field="ClassName",
                      out_table=r"/arcgis/home/output_table.gdb/output_table",
                      impervious_lyr.url)

Now we have the area of impervious surfaces per parcel, but only in a stand-alone table. Next, we'll join the stand-alone table to the Parcels attribute table. A table join updates the input table with the attributes from another table based on a common attribute field. Because we created the impervious area table with the name field from the Parcels layer, we'll perform the join based on that field.

Input
arcpy.management.JoinField(in_data=parcel_lyr.url, 
                           in_field="name", 
                           join_table="/arcgis/home/output_table.gdb/output_table", 
                           join_field="NAME")

Compute stormwater bill

We will now use add_to_definition method which adds a new field to an existing feature layer. We will then calculate percentage of impervious area to determine the tax value.

Input
parcel_lyr.manager.add_to_definition(json_dict={"fields":[{"name":"percent_imp",
                                                           "type":"esriFieldTypeDouble",
                                                           "alias":"percent_imp",
                                                           "nullable":True,
                                                           "editable":True}]})
parcel_lyr.manager.add_to_definition(json_dict={"fields":[{"name":"rate",
                                                           "type":"esriFieldTypeDouble",
                                                           "alias":"rate",
                                                           "nullable":True,
                                                           "editable":True}]})
Output
{'success': True}
Input
parcel_lyr.calculate(where="fid > 0", 
                     calc_expression=[{"field":"percent_imp",
                                       "sqlExpression":" ( impervious / ( pervious + impervious ) ) * 100"}])
parcel_lyr.calculate(where="fid > 0", 
                     calc_expression=[{"field":"rate",
                                       "sqlExpression":" impervious * 0.04"}])
Output
{'success': True,
 'updatedFeatureCount': 83142,
 'recordCount': 83142,
 'layerName': 'kent_county_parcels'}

percent_imp is the percentage of impervious area.

rate is the tax value calculated for $0.04 per square ft of impervious area.

Input
parcel_lyr.manager.delete_from_definition({"fields":[{"name":"name_1"}]})
Output
{'success': True}

Visualize results

Input
sdf = pd.DataFrame.spatial.from_layer(parcel_lyr)
sdf[['name', 'mailingadd', 'ownercity', 'propertyus', 'impervious', 'pervious', 'percent_imp', 'rate']]
Output
name mailingadd ownercity propertyus impervious_12 pervious_12 percent_imp rate
0 6-00-19700-01-1800-00001 OF HARRY RONALD WEBB GREENWOOD Agricultural-Land 330.0 645297.0 0.051113 13.20
1 6-00-19700-01-1802-00001 4318 MOUNT HOLLY RD GREENWOOD Vacant Land 0.0 47.0 0.000000 0.00
2 6-00-19700-01-1200-00001 25 ADAMSVILLE RD GREENWOOD Manuf Home/Ret.Title 2001.0 1773.0 53.020668 80.04
3 6-00-19600-01-5600-00001 143 LIDEN SCHOOL RD GREENWOOD Single-Fam Residential 935.0 6145.0 13.206215 37.40
4 6-00-19700-01-5900-00001 4299 MT HOLLY ROAD GREENWOOD Vacant Land 89.0 2226.0 3.844492 3.56
... ... ... ... ... ... ... ... ...
83137 1 291.0 94103.0 0.308282 11.64
83138 1 291.0 94103.0 0.308282 11.64
83139 1 291.0 94103.0 0.308282 11.64
83140 2-05-06818-01-2402-00001 2472.0 4624.0 34.836528 98.88
83141 2-05-06818-01-2400-00001 3500 S CLARK ST ARLINGTON Miscellaneous Imps. 19834.0 2270.0 89.730366 793.36

83142 rows × 8 columns

Symbolize the parcels

Input
m1 = gis.map('Deleware, kent county')
m1
Output

The parcels with the highest area of impervious surfaces are rendered with a darker color while the ones with less impervious surface are light in color. While we could symbolize the layer by the percentage of area that is impervious, most storm water fees are based on total area, not percentage of area. The code below applies this symbology using ClassedColorRenderer.

Input
m1.basemap = 'satellite'
Input
m1.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
                         'xmin': -8423844.205248041,
                         'ymin': 4736924.0084905885,
                         'xmax': -8386237.187331785,
                         'ymax': 4752211.414147602}
Input
m1.add_layer(lyr, { "type": "FeatureLayer",
                    "renderer":"ClassedColorRenderer",
                    "field_name":"impervious",
                    "opacity":0.7})

Conclusion

In this notebook, we reclassified an aerial image of a neighborhood in Louisville, Kentucky, to show areas that were pervious and impervious to water. We then determined the area of impervious surfaces per land parcel. With the information, the local government would be better equipped to determine storm water bills.

References

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.