A Methodology to Assess Land Use Development, Flooding, and Wetland Change as Indicators of Coastal Vulnerability

: Coastal areas around the world are becoming increasingly urban, which has increased stress to both natural and anthropogenic systems. In the United States, 52% of the population lives along the coast, and North Carolina is in the top 10 fastest growing states. Within North Carolina, the southeastern coast is the fastest growing region in the state. Therefore, this research has developed a methodology that investigates the complex relationship between urbanization, land cover change, and potential ﬂood risk and tested the approach in a rapidly urbanizing region. A variety of data, including satellite (PlanetScope) and airborne imagery (NAIP and Lidar) and vector data (C-CAP, FEMA ﬂoodplains, and building permits), were used to assess changes through space and time. The techniques consisted of (1) matrix change analysis, (2) a new approach to analyzing shorelines by computing adjacency statistics for changes in wetland and urban development, and (3) calculating risk using a ﬁshnet, or tessellation, where hexagons of equal size (15 ha) were ranked into high, medium, and low risk and comparing these results with the amount of urbanization. As other research has shown, there was a signiﬁcant relationship between residential development and wetland loss. Where urban development has yet to occur, most of the remaining area is at risk to ﬂooding. Importantly, the combined methods used in this study have identiﬁed at-risk areas and places where wetlands have migrated / transgressed in relationship to urban development. The combination of techniques developed here has resulted in data that local government planners are using to evaluate current development regulations and incorporating into the new long-range plan for the County that will include smart growth and identiﬁcation of risk. Additionally, results from this study area are being utilized in an application to the Federal Emergency Management Agency’s Community Response System which will provide residents with lower ﬂood insurance costs.


Introduction
Urbanization in the coastal area is an increasing worldwide problem where 10% of the total population and 13% of the urban population are located at 10 m or less above sea level [1]. In the United States, 52% of the population resides in coastal counties [2] and the Southeastern US has seen a rapid rate of population growth [3]. In the Southeastern US, North Carolina has experienced rapid population growth in the urban areas (i.e., Charlotte and Raleigh) and along the coast, where Brunswick, New Hanover, and Pender Counties have the fastest rates of population growth in the state [4]. Along with population growth, coastal flooding and storm-water management is an increasing problem in the Southeastern US including coastal North Carolina [5]. For example, Hurricanes Floyd (September 16,1999), Matthew (October 8, 2016), and Florence (September 12, 2018) resulted in catastrophic inland flooding to the region and very little flooding due to storm-surge.
Sea level has risen in North Carolina by 2.07-2.82 mm/year during the 20th century [6] and will likely rise 52-98 cm by 2100 [7], or 0.3 m at a low rate to an extreme rate of 25 m depending on the model used [8]. As sea level changes, the tidal plume may move further inland which leads to a change in wetland type. As such, when wetlands change, the ecosystem services that they provide will also change. The value of ecosystem services to the protection of urbanizing coastal areas has been well documented [9]. Theoretically, in flat terrain, salt marshes accrete sediment and transgress inland at the same rate as sea-level rise [10]. However, if sea-level rise surpasses the rate that marshes transgress inland, existing marshes will become intertidal mud flats or open water [11]. Further upstream, tidal freshwater wetlands will transition to salt tolerant vegetation as the plume moves further inland [12]. Freshwater wetlands accrete sediment at a faster rate than salt marsh vegetation. Thus, the potential transition from freshwater wetland to salt marsh vegetation would increase coastal susceptibility to inundation and flooding with increased sea-level rise predictions. In Southeastern North Carolina, planners recognize the importance of ecosystem services and the vulnerability of the coastal environment and need information to support responsible decision making to revise policies for managing urban development. Wetlands are sensitive and important natural resources that provide ecosystem services such as protection from storm impacts, but most research focuses on mapping wetlands and there is little information about the geographic connection between changes in wetlands and how these changes may relate to changes in upland areas that are becoming increasingly at risk to flooding [12][13][14]. Therefore, the purpose of this research was to utilize geospatial analysis techniques to synthesize a variety of data to derive indices of risk and to derive new measures of shoreline change with respect to wetland change and urban development. It was hypothesized that urban development has increased, wetland extent has decreased, and coastal flooding poses a risk to both urban and rural areas. Where wetlands have changed, we developed a method for classifying the type of change such as transgressing inland, becoming open water, or potential for flooding. Therefore, this research quantitatively investigated the spatial and temporal relationship between types of historic wetland change, flood potential, and urbanization in a manner that can be applied to other study areas. We utilized results to inform local government planning regulations and prepare for future land use development by identifying areas at risk to future development. Results from this research will help the local government to plan as they strive for sustainable development and coastal resilience. Additionally, methods utilized in this research may be translated to other coastal counties where flood potential, urban growth, and changes in wetlands are key factors for identifying coastal vulnerabilities for future urban development.

Study Area
Pender County is located in Southeastern North Carolina, on the South Atlantic coast of the United States ( Figure 1). The county is a mixture of agriculture, small towns, rivers, conservation lands, and a rapidly urbanizing coastal area along Highway 17. In general, Southeastern North Carolina has seen a substantial growth in population (Figure 2), where Brunswick and Pender counties had the fastest rates of growth (3.6% and 3.5%, respectively) from 2016 to 2017 [4].
A variety of data sources were gathered and analyzed (Table 1) and each data source is fully cited in the list of references. The data were analyzed using a variety of image processing and GIS techniques ( Figure 3). First, a digital elevation model (DEM) was assessed for potential stormwater flooding. Second, since Hurricane Florence made landfall in Wilmington, NC, on September 14, 2018 we analyzed the extent of flooding from this event. Third, land cover change analysis was conducted at two spatial scales (county-wide and coastal). Fourth, residential building through time was assessed using building permits. Lastly, results from these assessments were then integrated into a shoreline change analysis and urban development risk index.

Flooding and Potential Flooding
Using the 2014 North Carolina (NC) statewide QL2 Lidar data [22], a digital elevation model (DEM) for the Pender County study area (2,281 km 2 ) was generated using the ArcGIS 10.5.1 [23] terrain dataset tools according to the following steps: 1) convert LAS point cloud data [24] to a

Flooding and Potential Flooding
Using the 2014 North Carolina (NC) statewide QL2 Lidar data [22], a digital elevation model (DEM) for the Pender County study area (2,281 km 2 ) was generated using the ArcGIS 10.5.1 [22] terrain dataset tools according to the following steps: (1) convert LAS point cloud data [23] to a multipoint feature class, (2) create a feature dataset and import the LAS bare earth feature class, (3) build a terrain dataset, and (4) use the Terrain to Raster tool to generate the DEM [22]. The raw Lidar data was in the North Carolina State Plane coordinate system (FIPS 3200, NAD 83 datum and units in feet) and the resulting DEM had 5 × 5 ft cell size. The DEM was analyzed using flow direction and sink tools (ArcGIS 10.5.1) to detect topographic depressions where stormwater may collect. A sensitivity analysis was conducted using several kernel sizes, and it was decided that a minimum size of three adjacent cells (75 ft 2 ) removed smaller, and potentially erroneous, sinks while also keeping small areas that may collect stormwater.
Hurricane Florence made landfall in Wilmington, NC, on September 14, 2018 and subsequently produced record rainfall in Southeastern NC [24] during an already record-setting year for rainfall [25]. We used this opportunity to map the flooding from Hurricane Florence and compared this with depressions derived from the Lidar DEM to assess the validity of the geospatial model of potential flooding. PlanetScope 4-band (Red, Green, Blue and Near-Infrared) pre-and post-Hurricane Florence imagery were obtained from Planet [20]. Pre-Florence imagery was taken on August 15, 23, and 26, while Post-Florence imagery was taken on September 18 and 19. For both pre-and post-hurricane imagery, a normalized difference water index (NDWI) was calculated, and then an ISODATA unsupervised classification was conducted on the combined RGBN-IR and NDWI imagery. NDWI was calculated using this ratio: (Green-NearIR)/(Green+NearIR) [26].
Water was further processed to remove shadows (from trees and clouds) that were misclassified as water. The two classified images (pre and post storm) were compared to determine the extent of flooding attributable to the record rainfall from Hurricane Florence ( Figure 4). These maps were verified with local government employees who checked them against their records from the post hurricane damage assessment. unsupervised classification was conducted on the combined RGBN-IR and NDWI imagery. NDWI was calculated using this ratio: (Green-NearIR)/(Green+NearIR) [28].
Water was further processed to remove shadows (from trees and clouds) that were misclassified as water. The two classified images (pre and post storm) were compared to determine the extent of flooding attributable to the record rainfall from Hurricane Florence ( Figure 4). These maps were verified with local government employees who checked them against their records from the post hurricane damage assessment.

Land Cover Change
Two approaches were used to map land cover change through time. First, two land cover maps, from 1996 and 2010, from the National Oceanic Atmospheric Administration (NOAA)' Coastal Change Analysis Program (C-CAP) were used to identify land cover change [17]. The C-CAP program uses data from the National Landcover Dataset (NLCD) and has a spatial resolution of 30 ×

Land Cover Change
Two approaches were used to map land cover change through time. First, two land cover maps, from 1996 and 2010, from the National Oceanic Atmospheric Administration (NOAA)' Coastal Change Analysis Program (C-CAP) were used to identify land cover change [17]. The C-CAP program uses data from the National Landcover Dataset (NLCD) and has a spatial resolution of 30 × 30 m. Second, a higher spatial resolution approach was used in the coastal area to map wetlands from 1980 to 2016 using aerial photography (Aerial Photo Single Frames) and National Agriculture Imagery Program (NAIP) digital orthophoto quarter quads (DOQQs) [21,27]. Although some may consider aerial photography to be a dated source of information, the NAIP program provides useful imagery for mapping areas that are too large to be cost effective to purchase expensive high-resolution commercial imagery (e.g., WorldView), and this imagery also provides enough spectral information to derive general land cover classes. Therefore, while PlanetScope imagery worked well for county-wide analysis of flood inundation from Hurricane Florence, the PlanetScope imagery was tested with NAIP aerial photography for classifying land cover. NAIP provided more detail and better image classification than the PlanetScope imagery. Therefore, the NAIP imagery was used to map land cover along the coastal study area where greater detail was required ( Figure 5  The 2014 QL2 Lidar data were processed to create a normalized digital surface model (DSM) and DEM with 1 × 1 m spatial resolution to coincide with the NAIP 2016 imagery. The four bands (red (R), green (G), blue (B), and near-infrared (NIR)) from the NAIP imagery as well as the Lidar DSM and DEM layers were classified using eCognition [28]. For both image years (1980 and 2016), The 2014 QL2 Lidar data were processed to create a normalized digital surface model (DSM) and DEM with 1 × 1 m spatial resolution to coincide with the NAIP 2016 imagery. The four bands (red (R), green (G), blue (B), and near-infrared (NIR)) from the NAIP imagery as well as the Lidar DSM and DEM layers were classified using eCognition [28]. For both image years (1980 and 2016), an object-oriented classification technique using a multiresolution segmentation algorithm was used to segment the image into objects [29]. The objects were then classified using a decision tree process into five land cover types: water, emergent wetland, forest, agriculture, and developed land. Developed land was any area that has evidence of development (homes, roads, and managed grassland such as parks and planned open space). This decision-tree approach removed the speckled and potentially inaccurate representation of forested areas. An accuracy assessment was conducted where 50 randomly selected points, with a minimum of 30 m between points, were generated for each land cover class and each date/year. The overall accuracy for 1980 was 93% and 2016 was 92%. The two final maps where then compared using a classification change matrix. With two dates to compare, another approach to analyze change was a shoreline assessment of change, where each segment was assessed for the type of change that has occurred.

Urban Development
Residential building permits from 2006 through 2017 were obtained from Pender County [18] and were geocoded in ArcMap to analyze patterns of growth. The raw data were standardized to create addresses that were suitable for geocoding, and unmatched records were checked and re-matched until all acceptable data were geocoded. Some records did not have a street address and so these were not able to be geocoded; however, the iterative geocoding process resulted in an average of 90 percent of the records being geocoded (Table 2). To identify spatial patterns in the point data, a fishnet, or tessellation, using hexagons was generated for the coastal area. A variety of hexagon sizes were tested, and polygons of 15 hectares were chosen because this size was small enough to see patterns within subdivisions and not too small that there would be too many polygons with no data.

Coastal Vulnerability and Urban Development
To assess the relationship between urban development and vulnerability, several factors were aggregated to 15 ha tessellation polygons: total area of depressions, total area of wetlands, and average distance to floodplains. Each of these factors was converted to ranks (0, 1, 2, and 3) and then a total risk was calculated as the sum of the three ranked variables. Therefore, each tessellation polygon could have a total risk between 0 (all three variables ranked 0) and 9 (all three variables ranked with a value of 3). The total risk was then compared with the amount of developed land (measured as a percent of the tessellation polygon) to identify locations where there has been high urban development and high risk. Additionally, for areas that have not been fully developed (percent developed less than 25%), a comparison was made with risk to identify places that should be carefully considered before further development takes place. This last analysis is an important step in the spatial assessment of this coastal environment because it provides planners and decision-makers with a quantitative assessment of future development.

Flooding and Potential Flooding
The derived 2014 DEM for Pender County was analyzed to identify potential places where stormwater may collect. In total, 73 ha (97.44%) were identified as depressions within wetlands, forests, and grasslands while only 0.655 ha (0.09%) occurred within developed/urban areas (Table 3). Before Hurricane Florence, Pender County had 82.55 km 2 of surface water (mapped using August 2018 PlanetScope imagery), and a few days after the storm, an additional 123.27 km 2 of land were inundated (mapped using PlanetScope Imagery taken in September), resulting in 205.83 km 2 of inundated land (a 149% increase in area covered by water) (Figure 4). Flooding was greatest along the Cape Fear River, Northeast Cape Fear River, and Holly Shelter Game Land while there was much less flooding along the coast where there was an increase of only 0.17 ha (19%) ( Figure 6). Even though there was less flooding along the coast, this is where much of the population lives, which means the impact from flooding mostly impacted the rural part of the county, with possibly less access to emergency services, versus the coastal urban area which had less flooding. The potential flooding from the Lidar DEM analysis was compared with the actual flooding from Hurricane Florence, and although the total area is different (Hurricane Florence had much greater area of impact, 12,320 ha, than the small depressions generated using the DEM, 72.45 ha), there was no statistical difference between the locations of the two data products. For example, the areas of inundation from the hurricane and DEM depressions were intersected with the 2010 land cover data to compute percent inundation by land cover type. These percentages were compared using correlation (correlation coefficient = 0.98544), regression (R 2 = 0.9711), and F Test (P = 0.74286). These analyses show that the two datasets were very similar and were not statistically different.

Land Cover Change
From 1996 to 2010, Pender County had a 4.39% loss of forest, 1.02% loss of wetlands, and 6.57% growth in urban area along the coast (Figure 7 and Table 4). Additionally, while there was only a 1% loss in wetlands for all of Pender County, there was more (2.28%) loss of wetlands in the coastal zone where there was also high growth in urban development. Across Pender County, the loss of forest is attributable to an increase in agricultural/cleared land used for either timber or row crops (Table 5). Conversely, the higher resolution aerial photography (1980 and 2016) in the coastal area was dominated by 22.70% growth in urban development (from 12.69% to 35.40%) mostly because of a loss in forest (18.81%) and agriculture (6.86%) ( Table 6). Figure 8 illustrates the 2016 land cover (a) and the change to new urban land from 1980 to 2016 (b). Notice the large tracts of forest and agriculture that were converted to developed areas which are dominated by residential subdivisions. The potential flooding from the Lidar DEM analysis was compared with the actual flooding from Hurricane Florence, and although the total area is different (Hurricane Florence had much greater area of impact, 12,320 ha, than the small depressions generated using the DEM, 72.45 ha), there was no statistical difference between the locations of the two data products. For example, the areas of inundation from the hurricane and DEM depressions were intersected with the 2010 land cover data to compute percent inundation by land cover type. These percentages were compared using correlation (correlation coefficient = 0.98544), regression (R 2 = 0.9711), and F Test (P = 0.74286). These analyses show that the two datasets were very similar and were not statistically different.

Land Cover Change
From 1996 to 2010, Pender County had a 4.39% loss of forest, 1.02% loss of wetlands, and 6.57% growth in urban area along the coast (Figure 7 and Table 4). Additionally, while there was only a 1% loss in wetlands for all of Pender County, there was more (2.28%) loss of wetlands in the coastal zone where there was also high growth in urban development. Across Pender County, the loss of forest is attributable to an increase in agricultural/cleared land used for either timber or row crops (Table 5). Conversely, the higher resolution aerial photography (1980 and 2016) in the coastal area was dominated by 22.70% growth in urban development (from 12.69% to 35.40%) mostly because of a loss in forest (18.81%) and agriculture (6.86%) ( Table 6).      A statistical comparison between the land cover changes in the county versus the coast showed that there was low correlation (correlation coefficient = 0.49195) and low regression (R 2 = 0.3919), but an ANOVA determined that there was no statistical difference between the two change matrices. It is expected that the two change matrices (county and coast) should not be statistically different given their geographic proximity; however, there are clear differences between the two as described above. Therefore, an ANOVA comparison using only the changes (off-diagonals) showed that there was a A statistical comparison between the land cover changes in the county versus the coast showed that there was low correlation (correlation coefficient = 0.49195) and low regression (R 2 = 0.3919), but an ANOVA determined that there was no statistical difference between the two change matrices. It is expected that the two change matrices (county and coast) should not be statistically different given their geographic proximity; however, there are clear differences between the two as described above. Therefore, an ANOVA comparison using only the changes (off-diagonals) showed that there was a significant difference between the county and the coast (p = 0.078). The county was dominated by changes in forest while the coast had the largest change in urban growth, and even though the change in wetlands was small (1.02% in the county and 2.28% in the coast), the change in wetlands was significantly different along the coast in comparison with the inland part of the county. Regarding wetland change through time, wetlands comprised a smaller area compared to other cover types, which is expected. Hence, a closer look at the changes in coastal wetlands was necessary.
To quantify where wetland changes have taken place and how they have transitioned, we created a shoreline (total length of 143,774 m), divided it into equal lengths (1000 m), and classified each segment based on the type of wetland change along each segment (n = 145). Wetland change along the shoreline, from 1980 to 2016, had several types of change ( Figure 9):
Marsh migration inland/upland: wetland loss and adjacent wetland gain 3.
Sea-level rise and/or decrease in sedimentation: wetland loss with no change in urban areas 4.
Wetland gain and adjacent urban growth: increase in sedimentation resulted in the expansion of wetlands 5.
Wetland gain and no urban growth: increased sedimentation to maintain and expansion of wetlands Most of the study area shoreline (100 km or 69 percent) was class 1 (wetland loss and urban growth), which was expected given the large amount of urban development. Class 2 (wetland/salt marsh migration) was the next largest at 22 km (15 percent), and close behind was class 4 (wetland gain/urban growth) at 19 km (13 percent). Lastly, there was only 2 km (1.4 percent) of wetland gain (class 5) and 1 km (0.7 percent) of wetland loss with no other adjacent urbanization or wetland change (class 3). Urban growth: wetland loss and urban gain 2.
Marsh migration inland/upland: wetland loss and adjacent wetland gain 3.
Sea-level rise and/or decrease in sedimentation: wetland loss with no change in urban areas 4.
Wetland gain and adjacent urban growth: increase in sedimentation resulted in the expansion of wetlands 5.
Wetland gain and no urban growth: increased sedimentation to maintain and expansion of wetlands Most of the study area shoreline (100 km or 69 percent) was class 1 (wetland loss and urban growth), which was expected given the large amount of urban development. Class 2 (wetland/salt marsh migration) was the next largest at 22 km (15 percent), and close behind was class 4 (wetland gain/urban growth) at 19 km (13 percent). Lastly, there was only 2 km (1.4 percent) of wetland gain (class 5) and 1 km (0.7 percent) of wetland loss with no other adjacent urbanization or wetland change (class 3).

Urban Development
Geocoded new residential (not renovations or improvements) building permits from 2006 to 2017 were analyzed to identify patterns through space and time. After several iterations, 4537 records (90.31%) were geocoded ( Figure 10A). Both 2006 and 2007 had high levels of permits, and the county had more than the coast. However, these levels dropped considerably from 2008 through 2011, which mimics the 18 month national recession which began in December 2007. After the low period, building permits steadily increased from 2012 through 2015 and then declined through 2017. Unlike the earlier years, the coastal area had more permits from the rest of the county from 2012 to 2017.
The building permit points were summarized by polygon ( Figure 10B). The rate of change, or slope of the regression line, was calculated for each tessellation polygon which indicated growth (positive slope), no growth (slope near zero), or decline (negative slope) (Figure 11). Most hexagons (74.13%) did not have a positive or negative change in building permits from 2009 to 2017; however, all hexagons with the steepest positive slope were also located within the floodplain or within 2000 m of the floodplain, which indicates that the areas with the fastest growing urbanization also have the highest risk to flooding (

Urban Development
Geocoded new residential (not renovations or improvements) building permits from 2006 to 2017 were analyzed to identify patterns through space and time. After several iterations, 4537 records (90.31%) were geocoded ( Figure 10A). Both 2006 and 2007 had high levels of permits, and the county had more than the coast. However, these levels dropped considerably from 2008 through 2011, which mimics the 18 month national recession which began in December 2007. After the low period, building permits steadily increased from 2012 through 2015 and then declined through 2017. Unlike the earlier years, the coastal area had more permits from the rest of the county from 2012 to 2017.
The building permit points were summarized by polygon ( Figure 10B). The rate of change, or slope of the regression line, was calculated for each tessellation polygon which indicated growth (positive slope), no growth (slope near zero), or decline (negative slope) (Figure 11). Most hexagons (74.13%) did not have a positive or negative change in building permits from 2009 to 2017; however, all hexagons with the steepest positive slope were also located within the floodplain or within 2000 m of the floodplain, which indicates that the areas with the fastest growing urbanization also have the highest risk to flooding (Table 7).    * Distance of "0" indicates these building permits were located within the floodplain.
The land cover data (percent urban/developed in 2016 and the percent change in developed areas from 1980 to 2016, Figure 12) were compared with the total number of residential building permits and the rate/slope of residential building permits. Hypothetically, if there is a strong and direct relationship between the land cover data (derived from imagery) and the history of building permits, then if there is a shortage of imagery, geocoding building permits can be a useful secondary source of data to identify areas of growth. However, in this study area, there was no direct relationship between total number of permits and amount of development (R 2 = 0.1322) and no direct relationship between the rate of permits through time and the percent change in developed area (R 2 = 0.0404). Therefore, the two data sources are providing independent information about the urbanization of the coastal study area, where 1) the percentage of the tessellation polygon that is urban can relate to the degree of impervious surface and 2) the rate of change (slope) of permits indicates the direction of urbanization. Both are indicators of urbanization, but in this case, they are providing complimentary and independent information rather than redundant information.  * Distance of "0" indicates these building permits were located within the floodplain.
The land cover data (percent urban/developed in 2016 and the percent change in developed areas from 1980 to 2016, Figure 12) were compared with the total number of residential building permits and the rate/slope of residential building permits. Hypothetically, if there is a strong and direct relationship between the land cover data (derived from imagery) and the history of building permits, then if there is a shortage of imagery, geocoding building permits can be a useful secondary source of data to identify areas of growth. However, in this study area, there was no direct relationship between total number of permits and amount of development (R 2 = 0.1322) and no direct relationship between the rate of permits through time and the percent change in developed area (R 2 = 0.0404). Therefore, the two data sources are providing independent information about the urbanization of the coastal study area, where 1) the percentage of the tessellation polygon that is urban can relate to the degree of impervious surface and 2) the rate of change (slope) of permits indicates the direction of urbanization. Both are indicators of urbanization, but in this case, they are providing complimentary and independent information rather than redundant information.

Comparison of Urban Development with Flood Risk
Topographic depressions (identified from processing Lidar data), distance to floodplain, and total area of wetlands are the primary factors that indicate risk along the coast (Figure 13). Most of the coastal area has topographic depressions, which increases vulnerability to stormwater flooding, (Figure 13a). The distance to the 100 year floodplain was calculated (Figure 13b) and the total area of wetlands (Figure 13c) was the combination of total wetland area in the most recent imagery (2016) as well as the area of wetland loss (2008-2016) because risk is related to the existing wetland habitat as well as areas that are no longer wetlands since wetland habitats ameliorate the impact of storms. Total risk was computed by ranking (0, 1, 2, or 3) each of the three variables (area of depressions, distance to floodplain, and area of wetlands) and summing the ranks. Figure 13 illustrates the ranks and total risk. Of the total number of hexagons (n = 804), only three had 0 risk, while 318 (39.6%) had low risk (1, 2, or 3), 324 (40.3%) had medium risk (ranks 4 and 5), and 159 (20%) had high risk (ranks, 6, 7, 8, and 9). Only four hexagons had the highest possible ranking (9), which means these had the highest risk for all three factors.
Coastal resiliency necessitates an understanding of where urban development is occurring in relationship with the level of risk. Therefore, we compared risk with urbanization to identify where current development is at risk and where potential future development may be at risk and can be protected or land use development regulations can be adopted to protect homeowners and their property. The area of urban land in 2016 was summarized to the tessellation polygons, and although there is still undeveloped land, most of the coastal area (666 tessellations or 83%) has been at least 25% developed (Figure 14). These developing hexagons are a mixture of low risk (273 or 34%), medium risk (269 or 33%), and high risk (124 or 15%), which means 48% of the study area is urban and at risk. Conversely, of the hexagons that are not yet developed (urban ≤25%) (138 out of 804), there are 55 (7%) that have a medium risk and 35 (4%) that have a high risk, which means that 65% of the undeveloped area is at risk and may get developed, and these places need to have careful considerations to reduce the risk to these properties.

Comparison of Urban Development with Flood Risk
Topographic depressions (identified from processing Lidar data), distance to floodplain, and total area of wetlands are the primary factors that indicate risk along the coast (Figure 13). Most of the coastal area has topographic depressions, which increases vulnerability to stormwater flooding, (Figure 13a). The distance to the 100 year floodplain was calculated (Figure 13b) and the total area of wetlands (Figure 13c) was the combination of total wetland area in the most recent imagery (2016) as well as the area of wetland loss (2008-2016) because risk is related to the existing wetland habitat as well as areas that are no longer wetlands since wetland habitats ameliorate the impact of storms. Total risk was computed by ranking (0, 1, 2, or 3) each of the three variables (area of depressions, distance to floodplain, and area of wetlands) and summing the ranks. Figure 13 illustrates the ranks and total risk. Of the total number of hexagons (n = 804), only three had 0 risk, while 318 (39.6%) had low risk (1, 2, or 3), 324 (40.3%) had medium risk (ranks 4 and 5), and 159 (20%) had high risk (ranks, 6, 7, 8, and 9). Only four hexagons had the highest possible ranking (9), which means these had the highest risk for all three factors Coastal resiliency necessitates an understanding of where urban development is occurring in relationship with the level of risk. Therefore, we compared risk with urbanization to identify where current development is at risk and where potential future development may be at risk and can be protected or land use development regulations can be adopted to protect homeowners and their property. The area of urban land in 2016 was summarized to the tessellation polygons, and although there is still undeveloped land, most of the coastal area (666 tessellations or 83%) has been at least 25% developed ( Figure 14). These developing hexagons are a mixture of low risk (273 or 34%), medium risk (269 or 33%), and high risk (124 or 15%), which means 48% of the study area is urban and at risk. Conversely, of the hexagons that are not yet developed (urban ≤ 25%) (138 out of 804), there are 55 (7%) that have a medium risk and 35 (4%) that have a high risk, which means that 65% of the undeveloped area is at risk and may get developed, and these places need to have careful considerations to reduce the risk to these properties.

Discussion
The goal was to utilize a variety of readily available data and a combination of geospatial analysis techniques to investigate the relationship between land cover change, urbanization, wetland loss, and potential flooding with respect to how these factors spatially interact. Most studies of this type have focused on a relatively small study area, such as one urban area or a bay/estuary, where the investigation identifies where changes have occurred and computes various metrics (e.g., rate of Figure 14. Urbanization (percent) compared with risk: places with at least 25% development (2016) or increasing development (slope greater than 0.1) of medium risk (orange) or high risk (red) and places with low development (less than 25 percent) of medium (light green) or high risk (dark green).

Discussion
The goal was to utilize a variety of readily available data and a combination of geospatial analysis techniques to investigate the relationship between land cover change, urbanization, wetland loss, and potential flooding with respect to how these factors spatially interact. Most studies of this type have focused on a relatively small study area, such as one urban area or a bay/estuary, where the investigation identifies where changes have occurred and computes various metrics (e.g., rate of wetland loss) [30][31][32][33][34][35][36][37][38][39][40]. Conversely, in this project, we investigated a large study area (2281 km 2 ) where we conducted a generalized/regional assessment of the changing environment and investigated a large coastal area (105 km 2 ) in much greater detail using high-resolution data to look for patterns and trends for a variety of factors. Unlike previous studies, it is the combination of techniques that enabled the careful investigation of several changes throughout the area. We statistically tested the changes in the county versus the coast and identified that (1) the larger county area was dominated by a loss in forested area while the coast had the largest increase in urban development and (2) the change in wetlands was significantly greater along the coast in comparison with the inland part of the county.
Unlike previous research, this project utilized a vast amount of data at a very high resolution. For example, we analyzed aerial photography for the entire coastal zone including classifying the shoreline (144 km) into the type of changes that have taken place. We demonstrated that computing resources, software, data, and algorithms are fully capable of analyzing a large volume of data at multiple spatial scales to derive an understanding of the complex changes and the interplay among the various coastal processes at work in this area. Future research could build upon this multi-scale approach by using imagery obtained from unmanned aerial systems (UASs) at a very local scale where land cover and shoreline changes have been identified [41][42][43][44][45]. Instrumentation aboard unmanned aerial vehicles (UAVs) can provide excellent spatial resolution as well as multispectral bands; however, the coverage is limited because of battery time for overflights and access to launch and retrieve equipment can be problematic [44,45]. Therefore, a multi-scale approach that identifies critically changing locations is a prudent approach to identifying and ranking UAS missions.
The primary drivers of change to wetlands in a coastal environment are subsidence (which is minimal to non-existent in this study area), sea-level rise (SLR), and sedimentation. This project has documented that land cover change and urbanization along the coast has been substantial and is directly related to the natural and anthropogenic processes. Although the overall rate of change in wetlands was small (2.28% loss in wetlands along the coast), which is normal for a large study area, the location of change was significantly related to the processes that impact the coast. For example, most of the study area shoreline (100 km or 69%) had wetland loss and adjacent urban growth while a small amount (19 km or 13%) of shoreline had wetland gain in areas of urbanization which indicates that the sedimentation rate from development is keeping up with SLR. Evidence of wetland migration (22 km or 15%) demonstrated that SLR has necessitated that the environment keeps up with this rising water level, and subsequently, new wetlands/salt marsh have been established in the landward direction of the original marsh. These results confirm the conceptual model of wetland transition and migration [10,11] and the methods to derive this information. In this study area, urban expansion was the dominant detrimental impact to the coastal salt marsh; however, the variety of marsh changes provided examples of all types of coastal wetland change which demonstrates that the coastal environment is complex but directly correlated urbanization with wetland change.
The use of high spatial resolution Lidar data to derive a DEM, which was then used to identify potential places of storm-water collection/flooding, worked well in this study area. We utilized Planet satellite imagery to map water before and after a hurricane (Florence) to derive a map of floodwater, and when this was compared with DEM-derived potential stormwater collection points, these data matched well (significantly correlated) with the recent storm event. What this means is that there is value in having a historical record of flooding, but the additional topographic assessment based on a high spatial resolution DEM provides an independent predictive assessment of flood potential, which helps with planning and decision making to ameliorate stakeholders' concerns that may negate historical events. Therefore, although each site area is different and Lidar data vary considerably, it is a useful method to potentially identify stormwater management issues. Additionally, potential stormwater depressions provide an important addition to support local government planning and the management of rapidly urbanizing areas where stormwater management is an increasing problem in many areas [46][47][48][49][50][51]. Traditionally, FEMA floodplain boundaries are the de-facto standard to measure flood risk, but they have been highly criticized and do not provide a complete picture of flood potential [52,53].
The identification of coastal risk from flooding and wetland loss documented areas that can be further investigated. For example, future analysis could incorporate the Lidar elevation data (DSM and DEM) along with the wetland maps (which were mapped using NAIP aerial photography and Lidar) and change results to further identify areas at risk to wetland loss [54]. The most popular model to assess potential changes to coastal marshes due to SLR is SLAMM (sea-level affecting marshes model) [55]; however, there is an increasing amount of evidence that the output from this model is directly related to the accuracy (horizontal and vertical) of elevation and habitat mapping data [56]. Therefore, another potential future analysis with data derived from this project will be to test the accuracy of the SLAMM model in this study area. Now that we have documented areas at risk, it will be interesting and potentially useful to see if SLAMM also identifies these same areas at risk. Considering the reliance on elevation data (which can be inaccurate ), many researches have questioned the validity of the SLAMM output [54,56]. For example, one major drawback of SLAMM is that it does not incorporate urban development, or other land use conversions, into the model and since this is a driver of change in this study area, another modeling option is the use of ST-SIM (state-and-transition simulation model) [78]. The advantage of ST-SIM is that the model incorporates all types of changes, or transitions, and can simulate various scenarios such as SLR, urban development, loss of forest, agriculture, wetlands, etc. The simulations can then be compared to previous historical data to validate the transitions and plan for future land use development. We are currently testing ST-SIM at various locations along the North Carolina coast, such as barrier islands and large tidal creeks, where there have been changes in the tidal freshwater habitats that have been transitioning to saltwater marshes [12,79].
With regards to urban development, this study utilized an often unused and yet valuable data source: building permits. All property owners in the United States must purchase a building permit to develop land and therefore these records can be useful, if geocoded, to identifying spatial and temporal patterns. The study area described here has been rapidly urbanizing, like most coastal areas around the world; hence, most of the area is already developed or is in the process of been converted from natural to developed land. Only 11% of the study area is less than 25% developed and yet much of this area is at risk to possible flooding, wetlands loss, and short distance to floodplains. In this study, we used a tessellation approach to quantifying coastal risk and comparing it with urban development, but there are other approaches that can be utilized. For example, a multi-criteria stochastic modeling approach is simple and popular [80][81][82][83][84][85] but runs the risk of not accounting for inherent error in spatial data. Alternatively, fuzzy set methods and change vector analysis have demonstrated utility in the coastal zone [35,36,[86][87][88][89][90].

Conclusions
Analyzing coastal change is dependent on several factors: (1) natural and anthropogenic processes, (2) scale (spatial and temporal) of available data, and (3) techniques to identify and assess patterns and trends. In this study, we have conducted a multi-scale approach to identify coastal changes in a rapidly urbanizing area of southeastern North Carolina, USA. Urbanization is the largest driver of change where forest, agriculture, and wetlands have been replaced with urban development. Along the coast, 75% of the area and 69% of the shoreline has been urbanized while wetlands have been lost. Geocoded building permits and the calculation of change through time has identified areas of growth that provided new information not captured in the image processing and land cover change analysis.
Areas at risk were identified using Lidar data to derive depressions that are potential stormwater flooding, wetland loss which reduces the ability to ameliorate storm impacts, and distance to floodplains. These at-risk areas were compared with urbanization to reveal current problem areas and areas that are yet to be developed. These results are being used by planning officials to compare with existing planning regulations to propose new regulations to the Planning Board that will establish greater coastal resiliency through measures such as preserving wetlands, creating connected green spaces, and possibly adjusting building codes such as modifications to existing freeboard regulations which is more expensive to construct, but provides homeowners with a safer house if a flood event is to occur [91]. Additionally, Pender County is pursuing admission to FEMAs Community Rating System (CRS), which will enable county residents to participate in the National Flood Insurance Program (NFIP) at reduced rates [92,93]. However, to be granted into this program the county must meet planning standards that demonstrate that they are reducing risk to people and property. Information derived from this project is being used to apply to the CRS.
Lastly, this paper has demonstrated that the combination of image processing and geospatial analysis has identified patterns that would not be evident if only one approach had been undertaken. Investigating patterns using a variety of techniques will identify consistent patterns and specific locations where changes have occurred. Therefore, our recommendation is to utilize the data and methods described here to identify coastal changes and risk in other locations. Through this systematic approach to both raster and vector analysis, new information about coastal environments can be derived and compared. By utilizing multiple spatial analysis approaches, there can be verifiable measures of certainty and patterns in the landscape can be confirmed or reveal new information about how a place has been changing.