Spatial Assessment of Community Resilience from 2012 Hurricane Sandy Using Nighttime Light

Quantitative assessment of community resilience is a challenge due to the lack of empirical data about human dynamics in disasters. To fill the data gap, this study explores the utility of nighttime lights (NTL) remote sensing images in assessing community recovery and resilience in natural disasters. Specifically, this study utilized the newly-released NASA moonlight-adjusted SNPP-VIIRS daily images to analyze spatiotemporal changes of NTL radiance in Hurricane Sandy (2012). Based on the conceptual framework of recovery trajectory, NTL disturbance and recovery during the hurricane were calculated at different spatial units and analyzed using spatial analysis tools. Regression analysis was applied to explore relations between the observed NTL changes and explanatory variables, such as wind speed, housing damage, land cover, and Twitter keywords. The result indicates potential factors of NTL changes and urban-rural disparities of disaster impacts and recovery. This study shows that NTL remote sensing images are a low-cost instrument to collect near-real-time, large-scale, and high-resolution human dynamics data in disasters, which provide a novel insight into community recovery and resilience. The uncovered spatial disparities of community recovery help improve disaster awareness and preparation of local communities and promote resilience against future disasters. The systematical documentation of the analysis workflow provides a reference for future research in the application of SNPP-VIIRS daily images.


Introduction
Natural disasters continue to cause widespread and long-lasting economic, social, and environmental impacts during the past decades. Hurricanes, one of the devastating natural hazards in the United States, adversely interrupt businesses, impact transportation, and disrupt communities. Empirical evidence shows that communities may suffer from different consequences and recover at different rates during a disaster. Such differences can be attributed to community resilience, which is the critical ability of individuals and communities to cope with, adapt to, and recover from external stresses [1]. To reduce the adverse impacts of disasters, considerable efforts have been made to understand and measure community resilience. One of the challenges in this endeavor is the lack of empirical data to monitor different aspects of human activities in disasters [2]. Previous studies on community resilience mostly rely on qualitative data (e.g., surveys and interviews) [3][4][5] and resilience indices aggregated from socio-economic variables [6][7][8][9][10]. However, the qualitative data lack timeliness, are costly to acquire, and are inapplicable in large geographic areas [2,11]. The resilience index-based approaches were theory-driven and lack validation from empirical data [2,12]. Thus, there is a pressing need for large-scale, high-resolution data to detect human dynamics in disasters, from which to extract new indicators to measure community resilience.
With the development of remote sensing and earth observation techniques, various types of images of the earth's surface are captured at different spatial and temporal resolutions. As a special type of remote sensing techniques, nighttime light (NTL) remote sensing Hurricane Sandy (2012), also called Superstorm Sandy, was formed in the western Caribbean Sea on 22 October, made its landfall near Atlantic City in New Jersey in the United States on 29 October, and finally dissipated on 2 November. Hurricane Sandy impacted the Northeastern United States with heavy rain and strong winds and later caused flash flooding and storm surge consequently. In total, 233 people were dead across eight countries from the Caribbean to Canada and over 650,000 homes were damaged in the US [34]. Moreover, Hurricane Sandy caused around 68.7 billion worth of damage, which heavily impacted the economy in the Northeastern US. According to the Hazus disaster damage model from the Federal Emergency Management Agency (FEMA), 12 states were heavily impacted by Hurricane Sandy, including New York, Connecticut, New Jersey, Delaware, District of Columbia, Maryland, Massachusetts, New Hampshire, Pennsylvania, Rhode Island, Virginia, and West Virginia [35]. To capture human dynamics in populated areas, 66 Core-Based Statistical Areas (CBSAs) within 12 states, including both metropolitan and micropolitan statistical areas, are selected as the study area (yellow areas in Figure 1) (CBSA data can be downloaded from the US Census [36]). According to 2010 American Community Survey (2012ACS 5-year estimates) and Bureau of Economic Analysis (BEA) data, the study area has a 19.7% of the US total population (60,876,820 out of 308,745,538) and contributes 28.7% of the Gross Domestic Product generated in the US in 2012 ($4,107,550,656,000 out of $14,332,171,020,000). Within these 66 CBSAs, 52.3 million people (86.0% of the total population) lived in coastal areas. The population in the study area consists of 69.0% White, 16.1% African American, and 6.2% Asian. The research area has 5 of the 25 largest urbanized areas in the US, including Washington D.C., Philadelphia, New York City, Boston, and Baltimore [37]. In summary, the study area is one of the most urbanized regions in the US and it is of utmost importance to analyze how a well-developed region responds and recovers from one of the costliest Atlantic hurricanes in US history. Four SNPP-VIIRS DNB tiles (h09v05, h10v04, h10v05, and h11v04) that fully cover the study area were utilized to analyze NTL radiance changes. The spatial analysis of NTL radiance has been conducted in pixels, block groups, and CBSAs to avoid the Modifiable Areal Unit Problem (MAUP).

NASA Black Marble Products
As the primary data source in this study, NASA Black Marble Product, also known as the VNP46 product suite, includes two types of images: (1) NPP/VIIRS At-sensor TOA Nighttime Radiance Daily L3 Global 500 m SIN Grid V001 (NASA product ID: VNP46A1)

NASA Black Marble Products
As the primary data source in this study, NASA Black Marble Product, also known as the VNP46 product suite, includes two types of images: (1) NPP/VIIRS At-sensor TOA Nighttime Radiance Daily L3 Global 500 m SIN Grid V001 (NASA product ID: VNP46A1) and (2) NPP/VIIRS Moonlight-adjusted Nighttime Lights Daily L3 Global 500 m SIN Grid V001 (NASA product ID: VNP46A2). VNP46A2 products are enhanced from VNP46A1 by screening moonlight and calibrating the BRDF. VNP46A1 products only contain fundamental layers, such as spectral bands, acquisition time, and original radiance collected from the sensor. Improved from VNP46A1, the VNP46A2 products include additional layers, such as BRDF-corrected radiance, BRDF-corrected gap-filled radiance, mandatory quality assurance (QA) flags, and snow flags, which inform the data quality and correction been conducted at each pixel. Images used in this study were downloaded from NASA's Level 1 and Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS-DAAC) using the modified PowerShell wget tool (version 5.1.18362.1801) in Windows OS with two parameters: (1) Tiles: four tiles (h09v05, h10v04, h10v05, and h11v04) covering the study area (Figures 1 and 2 (value 00) in the non-filled images (VNP46A2) and excluded pixels in the other quality levels (01, 02, and 255). Daily NTL images of the research area that are collected by a single scan and mainly covered by high-quality pixels were used. Then the qualified daily NTL images were mosaiced in R and later cropped into different spatial units (e.g., block groups, CBSAs, etc.) for spatial and explanatory analysis.

Auxiliary Data
Auxiliary data were collected from multiple sources to explain the variation of NTL during Hurricane Sandy (Table 1). The auxiliary data include wind speed data, damaged housing unit data, Twitter data, land-use data, and proximity data, all of which are considered to have potential relations with NTL variation. Wind speed data at ground observations in the study area were obtained from the National Oceanic and Atmospheric Administration (NOAA) Hurricane Sandy aftermath report [39]. Damaged housing unit data at the block-group level were acquired from the U.S. Department of Housing and Urban Development (HUD) [40]. Twitter data were collected from Archive.org (an online archive of digital materials) from 22 October 2012 to 17 November 2012 and aggregated at the Core-Based Statistical Area (CBSA) scale according to geotags and user locations. Defined by the Office of Management and Budget (OMB), a CBSA is a geographic area that consists of one or more counties (or equivalents) anchored by an urban center of at least 10,000 people plus adjacent counties that are socioeconomically tied to the urban center by commuting [41]. CBSA refers collectively to both Metropolitan Statistical Areas (MSA) and Micropolitan Statistical Areas (μSA). The NLCD 2011 Land Cover (product name: NLCD_2011_Land_Cover_L48_20190424) and NLCD 2011 Developed Imperviousness Descriptor (product name: NLCD_2011_Impervious_descriptor_L48_20190405) in the Contiguous United States (CONUS) from the Multi-Resolution Land Characteristics (MRLC) Consortium (http://www.mrlc.gov) accessed on 18 September 2021 [42][43][44][45][46][47] were used to compare NTL changes in different land uses. Spatial analysis tools were applied to aggregate and analyze the NTL data in different spatial units and land-use types. In addition, the hurricane trajectory was collected from National Weather Service (NWS) and was used to calculate the proximity to Hurricane Sandy using Euclidean Distance. The NTL radiance data were extracted from the compressed Black Marble images in Hierarchical Data Format (HDF5) using R. The BRDF-corrected radiance (in VNP46A1 products), BRDF-corrected gap-filled radiance (VNP46A2), and non-corrected radiance at sensor (in both VNP46A1 and VNP46A2) were compared in Figure 2. The result shows that the corrected radiance (VNP46A2 & VNP46A2 gap-filled) has a more stable daily NTL radiance compared with the original radiance (VNP46A1) that is fluctuated due to moonlight cycles ( Figure 2). The UTC_Time layer in the VNP46A1 products contains the UTC time when nighttime radiance is acquired. Four types of QA values in the VNP46A2 products indicate the quality of pixels: value 00 denotes high-quality persistent NTL, value 01 means high-quality ephemeral NTL, value 02 is poor quality (outlier, potential cloud contamination, or other issues), and value 255 is no retrieval (VNP46A2 gap-filled layer applies high-quality pixels from the nearest date to replace no-retrieval pixels) [32,38]. To ensure high-quality and consistent NTL radiance, this study only used high-quality pixels (value 00) in the non-filled images (VNP46A2) and excluded pixels in the other quality levels (01, 02, and 255). Daily NTL images of the research area that are collected by a single scan and mainly covered by high-quality pixels were used. Then the qualified daily NTL images were mosaiced in R and later cropped into different spatial units (e.g., block groups, CBSAs, etc.) for spatial and explanatory analysis.

Auxiliary Data
Auxiliary data were collected from multiple sources to explain the variation of NTL during Hurricane Sandy ( Table 1). The auxiliary data include wind speed data, damaged housing unit data, Twitter data, land-use data, and proximity data, all of which are considered to have potential relations with NTL variation. Wind speed data at ground observations in the study area were obtained from the National Oceanic and Atmospheric Administration (NOAA) Hurricane Sandy aftermath report [39]. Damaged housing unit data at the block-group level were acquired from the U.S. Department of Housing and Urban Development (HUD) [40]. Twitter data were collected from Archive.org (an online archive of digital materials) from 22 October 2012 to 17 November 2012 and aggregated at the Core-Based Statistical Area (CBSA) scale according to geotags and user locations. Defined by the Office of Management and Budget (OMB), a CBSA is a geographic area that consists of one or more counties (or equivalents) anchored by an urban center of at least 10,000 people plus adjacent counties that are socioeconomically tied to the urban center by commuting [41]. CBSA refers collectively to both Metropolitan Statistical Areas (MSA) and Micropolitan Statistical Areas (µSA). The NLCD 2011 Land Cover (product name: NLCD_2011_Land_Cover_L48_20190424) and NLCD 2011 Developed Imperviousness Descriptor (product name: NLCD_2011_Impervious_descriptor_L48_20190405) in the Contiguous United States (CONUS) from the Multi-Resolution Land Characteristics (MRLC) Consortium (http://www.mrlc.gov, accessed on 18 September 2021) [42][43][44][45][46][47] were used to compare NTL changes in different land uses. Spatial analysis tools were applied to aggregate and analyze the NTL data in different spatial units and land-use types. In addition, the hurricane trajectory was collected from National Weather Service (NWS) and was used to calculate the proximity to Hurricane Sandy using Euclidean Distance.

Conceptual Framework of Recovery Trajectories
Based on the concept of resilience including resistance and recovery, resilience under the disastrous condition should be measured on both the reduction of functional capacities and the recovery of the capacities to a normal condition [48]. Recovery trajectories conceptualize resilience as a dynamic process that describes the change of a social system's functional capacities after a shock, such as natural disasters. Figure 3 shows an example of the recovery trajectory where a functional capacity suddenly declines after the disaster and gradually recovers subsequently. The conceptual framework of recovery trajectories has been widely used to measure the resilience of social-ecological and infrastructural systems [49][50][51][52][53]. The variation of recovery trajectories among different places can indicate resilience. In this study, due to varying the atmospheric condition during the hurricane, many images do not have sufficient high-quality pixels covering the study area. The time series of the images cannot provide continuous recovery trajectories. Thus, the analysis of NTL changes is limited to three sampling points at pre-disaster, in-disaster, and post-Remote Sens. 2021, 13, 4128 6 of 20 disaster phases. Specifically, the Disturbance Rate (D NTL ) and Recovery Rate (R NTL ) of NTL radiance are measured from the three sampling points: where NTL pre , NTL in , and NTL post denote the NTL radiance measured in the pre-, in-, and post-disaster phases respectively ( Figure 3). The function E calculates the mean NTL radiance within a spatial unit (e.g., block groups or CBSAs) or land use type. In general, a small disturbance (low D NTL ) and fast recovery (high R NTL ) would indicate high resilience, while high D NTL and low R NTL mean the opposite.

NTL Image Selection and Processing
NTL data from four tiles were extracted from the VNP46 suite and mosaiced, including the UTC time layer, mandatory QA flag layer, and three radiance layers. Boundaries of 66 CBSAs in the Northeastern US were used as the extent of the extracted NTL data. Due to varying cloud cover, the data quality of NTL images is very unstable during the disaster period. Thus, we selected images with sufficient high-quality pixels in three disaster phases (pre-disaster, in-disaster, and post-disaster) to analyze NTL changes.
Referring to previous studies [33,54], we split the hurricane process into three periods by the landing time (29 October) and dissipating time (2 November), including the predisaster period (22)(23)(24)(25)(26)(27)(28), the in-disaster period (29 October-1 November), and the post-disaster period (2-17 November). Three sampling time points were used to represent NTL radiance in the three periods respectively. The selection is based on two criteria: data quality and single acquisition time. For the former criterion, dates are selected when the study area is mostly covered by high-quality (cloudless) pixels. For the latter criteria, only days when the study area is scanned at the same time are selected. Figure 4 shows the UTC_Time layers (denoting the acquisition time in UTC) in the VNP46A1 products from 22 October to 17 November in 2012. The blue dashed lines indicate days when the study area is covered by image tiles, scanned at different time points. As an example, on 8 November the study area is covered by two image tiles scanned at different time points (red and cyan in Figure 5a). A clear radiance difference can be observed in the boundary between the two tiles (see Figure 5b). Thus, mosaicking of images with various scanning time points on the same day can cause inconsistent NTL radiance over space. To avoid this bias, dates with multiple acquisition times (highlighted in blue dashed lines in Figure  6) were excluded from the selection.

NTL Image Selection and Processing
NTL data from four tiles were extracted from the VNP46 suite and mosaiced, including the UTC time layer, mandatory QA flag layer, and three radiance layers. Boundaries of 66 CBSAs in the Northeastern US were used as the extent of the extracted NTL data. Due to varying cloud cover, the data quality of NTL images is very unstable during the disaster period. Thus, we selected images with sufficient high-quality pixels in three disaster phases (pre-disaster, in-disaster, and post-disaster) to analyze NTL changes.
Referring to previous studies [33,54], we split the hurricane process into three periods by the landing time (29 October) and dissipating time (2 November), including the predisaster period (22)(23)(24)(25)(26)(27)(28), the in-disaster period (29 October-1 November), and the post-disaster period (2-17 November). Three sampling time points were used to represent NTL radiance in the three periods respectively. The selection is based on two criteria: data quality and single acquisition time. For the former criterion, dates are selected when the study area is mostly covered by high-quality (cloudless) pixels. For the latter criteria, only days when the study area is scanned at the same time are selected. Figure 4 shows the UTC_Time layers (denoting the acquisition time in UTC) in the VNP46A1 products from 22 October to 17 November in 2012. The blue dashed lines indicate days when the study area is covered by image tiles, scanned at different time points. As an example, on 8 November the study area is covered by two image tiles scanned at different time points (red and cyan in Figure 5a). A clear radiance difference can be observed in the boundary between the two tiles (see Figure 5b). Thus, mosaicking of images with various scanning time points on the same day can cause inconsistent NTL radiance over space. To avoid this bias, dates with multiple acquisition times (highlighted in blue dashed lines in Figure 6) were excluded from the selection.     Combining the two criteria, 22 October was selected as the sampling time point in the pre-disaster period. 17 November was selected as the sampling point in the post-disaster period with additional considerations. First, dates between 6 November and 10 November were excluded from the selection to eliminate impacts from the 2012 Nor'easter Storm. 12 November meets both criteria but does not cover the populated areas (e.g., New York City) near the landfall location, and thus was not selected. Due to the thick cloud brought by the hurricane, it has been a challenge to find a single image meeting both criteria during the in-disaster period. Thus, the two images on 31 October and 1 November were combined to create a larger spatial coverage. When combining the images, if a pixel had a high-quality flag on both days, the average NTL radiance of the two days was assigned in the combined image. If a pixel had a good QA value in one of the two days, the NTL radiance on that day was selected. Pixels with poor QA values on both days were excluded from the analysis.

Spatial Analysis
The spatial variation of NTL changes metrics ( and ) in Hurricane Sandy was analyzed in different spatial units, including the pixel, block group, and CBSA level, thus avoiding the MAUP in the spatial analysis. At the pixel level, and were calculated in high-quality pixels at the three sampling time points using Equations (1) and (2). For visualization purposes, a 5-by-5 focal window was applied to smooth the spatial pattern at the pixel level. At the block-group and CBSA level, zonal statistic methods were used to aggregate and in each spatial unit. To mainly focus on the illuminated areas at night, a mask of impervious surface area (NLCD 2011 Developed Imperviousness Descriptor) was applied to the NTL radiance before aggregating and . Then, choropleth maps were produced to visualize the spatial distribution of and at two levels. Comparing different distribution maps, spatial variations of NTL change as well as the underlying factors were discussed. Additionally, the zonal statistics tool was applied to extract and compare NTL radiance in different land uses. Combining the two criteria, 22 October was selected as the sampling time point in the pre-disaster period. 17 November was selected as the sampling point in the post-disaster period with additional considerations. First, dates between 6 November and 10 November were excluded from the selection to eliminate impacts from the 2012 Nor'easter Storm. 12 November meets both criteria but does not cover the populated areas (e.g., New York City) near the landfall location, and thus was not selected. Due to the thick cloud brought by the hurricane, it has been a challenge to find a single image meeting both criteria during the in-disaster period. Thus, the two images on 31 October and 1 November were combined to create a larger spatial coverage. When combining the images, if a pixel had a high-quality flag on both days, the average NTL radiance of the two days was assigned in the combined image. If a pixel had a good QA value in one of the two days, the NTL radiance on that day was selected. Pixels with poor QA values on both days were excluded from the analysis.

Spatial Analysis
The spatial variation of NTL changes metrics (D NTL and R NTL ) in Hurricane Sandy was analyzed in different spatial units, including the pixel, block group, and CBSA level, thus avoiding the MAUP in the spatial analysis. At the pixel level, D NTL and R NTL were calculated in high-quality pixels at the three sampling time points using Equations (1) and (2). For visualization purposes, a 5-by-5 focal window was applied to smooth the spatial pattern at the pixel level. At the block-group and CBSA level, zonal statistic methods were used to aggregate D NTL and R NTL in each spatial unit. To mainly focus on the illuminated areas at night, a mask of impervious surface area (NLCD 2011 Developed Imperviousness Descriptor) was applied to the NTL radiance before aggregating D NTL and R NTL . Then, choropleth maps were produced to visualize the spatial distribution of D NTL and R NTL at two levels. Comparing different distribution maps, spatial variations of NTL change as well as the underlying factors were discussed. Additionally, the zonal statistics tool was applied to extract and compare NTL radiance in different land uses.

Regression Analysis
Univariate ordinary least square regression was applied to examine the relation between NTL radiance changes (D NTL and R NTL ) and each explanatory variable. The equations of the regression model are as follows: where x i is the explanatory variable, y is D NTL or R NTL derived from NTL changes using Equations (1) and (2), and ε is the residual. First, wind speed data were analyzed at the 1-km buffer scale. The maximum surface wind speed data during Hurricane Sandy in the research area were collected from the NOAA Hurricane Sandy aftermath report, geocoded based on the geographic coordinates, and later used to explain the spatial variation of disturbance (D NTL ) and recovery of NTL (R NTL ) [39]. In total, 387 observation stations located in high-quality pixels were used to analyze the relations. Maximum surface wind speed occurred near the landfall time (from 29 October, 12:51 pm UTC to 30 October, 4:57 am UTC) at these stations were compared with NTL changes in 1-km buffers around the stations.
Second, damaged housing unit data were analyzed at the block-group scale. Data collected from HUD include block groups in New York, New Jersey, Connecticut, and Rhode Island. Percentage of non-seasonal damaged housing units in block groups were used (pct_dmg_ns in the HUD dataset) as an explanatory variable. This percentage was calculated using the number of damaged non-seasonal housing units divided by the number of total non-seasonal housing units.
Third, geocoded Twitter data were assessed at the CBSA scale. According to user profiles and geotags, the collected tweets were geocoded at the city level and linked to each affiliated CBSA. Hurricane-related tweets were retrieved by querying the keywords 'hurricane' and 'sandy'. The retrieved tweets were preprocessed by removing stop words (e.g., a, the, is, etc.) and word stemming (e.g., converting flooded/flooding/floods to flood). The average sentiments of hurricane-related and all tweets in each CBSA were calculated using the Valence Aware Dictionary and sEntiment Reasoner (VADER) package (version 3.2.1) in Python 3.8.2 (https://pypi.org/project/vaderSentiment/, accessed on 18 September 2021) [55,56]. The ratio of hurricane-related tweets to all tweets was calculated in CBSAs. Next, tweets containing various keywords (e.g., 'electric', 'close', 'outage', 'flood', 'evacuate', 'wind', 'rain', 'storm', and 'damage') were extracted from the hurricanerelated tweets. The selection of keywords refers to previous studies that utilize Twitter data to study natural disasters [54,57,58]. The ratio of tweets containing each keyword was calculated in CBSAs, with the number of tweets containing the keyword as the nominator and the number of all geocoded tweets as the denominator. The relation between tweets and NTL change was compared among CBSAs using univariate regression models.
Fourth, Euclidean Distance was applied to calculate the distance to the hurricane trajectory in each pixel. The relations of D NTL and R NTL in each pixel to distance to hurricane track were examined using regression analysis. Additionally, D NTL and R NTL . were compared in different land use and land cover (LULC) types. The comparison was conducted in two LULC classifications, including NLCD 2011 Land Cover (CONUS) and NLCD 2011 Developed Imperviousness Descriptor (CONUS). Zonal statistics were used to calculate the mean value of D NTL and R NTL in each land-use type.
Finally, results from all univariate regression models at different research scales are summarized in the table. The coefficients of these models (e.g., p-value, R 2 ) are reported to indicate the significance of the relationships. Variables with high correlation with D NTL and R NTL are highlighted. The potential driving factors of these significant correlations with the NTL change are discussed. The workflow of the entire analysis process is summarized in Figure 7.

Temporal Changes of NTL Radiance
Within the study area, 39 CBSAs with sufficient high-quality pixels were plotted in

Temporal Changes of NTL Radiance
Within the study area, 39 CBSAs with sufficient high-quality pixels were plotted in

Spatial Patterns of NTL Change
The spatial patterns of D NTL and R NTL are shown at the pixel and CBSA levels in Figure 9. At the pixel level (Figure 9a,b), the New York/New Jersey Bight area has a large NTL disturbance (high absolute value of D NTL , red color) in the hurricane and a strong recovery after the hurricane (high R NTL ). This area is highly populated and close to the landfall location of Hurricane Sandy. As an exception, the center of New York City only had moderate NTL disturbance and recovery. Additionally, NTL disturbance can also be observed in Providence, Virginia Beach. At the CBSA level (Figure 9c,d), the metropolitan areas near New York City and Providence had a large NTL reduction (high absolute D NTL ), but the NTL in most areas soon recovered from the reduction (high R NTL ). Regions near Providence, Philadelphia, New York, Washington D.C., have recovered to the pre-disaster condition and can be considered resilient from this perspective. Several inland areas, such as Lancaster and York-Hannover CBSA in Pennsylvania, experienced NTL reduction and did not recover to their pre-disaster condition. As can be observed in Figure 9, the NTL disturbance shows different spatial patterns when aggregated in different spatial units. For example, areas near Virginia Beach had a large NTL reduction at the pixel level but had a small reduction at CBSA level. The variation in different spatial levels demonstrates the MAUP when analyzing the NTL changes.

Spatial Patterns of NTL Change
The spatial patterns of and are shown at the pixel and CBSA levels in Figure 9. At the pixel level (Figure 9a,b), the New York/New Jersey Bight area has a large NTL disturbance (high absolute value of , red color) in the hurricane and a strong recovery after the hurricane (high ). This area is highly populated and close to the landfall location of Hurricane Sandy. As an exception, the center of New York City only had moderate NTL disturbance and recovery. Additionally, NTL disturbance can also be observed in Providence, Virginia Beach. At the CBSA level (Figure 9c,d), the metropolitan areas near New York City and Providence had a large NTL reduction (high absolute ), but the NTL in most areas soon recovered from the reduction (high ). Regions disaster condition and can be considered resilient from this perspective. Several inland areas, such as Lancaster and York-Hannover CBSA in Pennsylvania, experienced NTL reduction and did not recover to their pre-disaster condition. As can be observed in Figure  9, the NTL disturbance shows different spatial patterns when aggregated in different spatial units. For example, areas near Virginia Beach had a large NTL reduction at the pixel level but had a small reduction at CBSA level. The variation in different spatial levels demonstrates the MAUP when analyzing the NTL changes.

Univariate Regression Results
Univariate regression analyses were conducted to associate the NTL variation with explanatory variables. The purpose of the regression analysis is to find underlying factors that can explain the spatial variation of NTL radiance observed in the spatial analysis (Section 4.2). Depending on the resolutions of the auxiliary data, regression results at all spatial scales are listed in Table 2. In this study, we consider p < 0.05 and R 2 > 0.13 [59] as significant relations. Highlighted in bold font in Table 2, 6 of the total 36 regression models meet both criteria. First, a significant relation between NTL fluctuation ( and ) and wind speed is detected, indicating high wind speed caused a V-shape recovery pattern (low negative and high positive ) during the hurricane. This result

Univariate Regression Results
Univariate regression analyses were conducted to associate the NTL variation with explanatory variables. The purpose of the regression analysis is to find underlying factors that can explain the spatial variation of NTL radiance observed in the spatial analysis (Section 4.2). Depending on the resolutions of the auxiliary data, regression results at all spatial scales are listed in Table 2. In this study, we consider p < 0.05 and R 2 > 0.13 [59] as significant relations. Highlighted in bold font in Table 2, 6 of the total 36 regression models meet both criteria. First, a significant relation between NTL fluctuation (D NTL and R NTL ) and wind speed is detected, indicating high wind speed caused a V-shape recovery pattern (low negative D NTL and high positive R NTL ) during the hurricane. This result suggests that large D NTL disturbance occurred in areas that experienced severe physical impacts from the hurricane. At the CBSA scale, regression analysis was conducted to associate NTL changes with Twitter data. Significant relations are detected between: (1) D NTL and the ratio of keyword 'damage', (2) D NTL and the ratio of 'sandy', (3) R NTL and the ratio of 'electric'. The first relation suggests that the NTL disturbance is a signal of property damage. The significant relation with the keyword "sandy" indicates areas experienced high NTL disturbance also have more people discussing the hurricane in Twitter. The relation between R NTL and keyword 'electric' may imply that electric outage or restoration is an important factor in the recovery process. These detected relations indicate that the NTL changes may reflect certain aspects of human activities related to the Twitter keywords. Moreover, a significant relation between D NTL and distance to the hurricane track was detected at the CBSA level, indicating that a higher NTL reduction was observed near the hurricane track. Due to low p-value or R 2 , the other 30 regression models do not show significant relations.
Finally, NTL metrics (D NTL and R NTL ) were compared among different LULC types. Figure 10a shows that non-urban land cover types (open water, barren lands, cultivated crops) have a large radiance reduction (large absolute value of D NTL ) and relatively slow recovery (lower R NTL ). Urban areas, including high-, medium-, and low-intensity developed land, and open space had relatively small radiance reduction (small absolute D NTL ) but faster recovery (higher R NTL ). Among the four urban land uses, high-intensity developed area had the largest NTL reduction (largest absolute D NTL ) in the hurricane. In different imperviousness surfaces, primary and secondary roads had large radiance reduction (large absolute D NTL ) which may be caused by interrupted traffic (Figure 10b). The recovery rate of primary roads in urban areas is higher than primary roads in non-urban areas. Tertiary and secondary roads outside urban areas had little fluctuations (low absolute D NTL and low R NTL ). In general, roads within urban areas have shown a strong post-disaster recovery, while roads in rural areas take longer to recover.

Discussion
The analysis aims to evaluate the utility of NTL images as a potential data source to study community resilience. Compared with previous studies, this study has the following contributions. First, it demonstrates the utility of the new moonlight-adjusted Black Marble images (VNP46A2) in detecting short-term human dynamics in a disaster. Despite the advancements claimed by NASA, the utility of the new image product has not been

Discussion
The analysis aims to evaluate the utility of NTL images as a potential data source to study community resilience. Compared with previous studies, this study has the following contributions. First, it demonstrates the utility of the new moonlight-adjusted Black Marble images (VNP46A2) in detecting short-term human dynamics in a disaster. Despite the advancements claimed by NASA, the utility of the new image product has not been tested in real disaster events. This study systematically documented the analytical workflow and noted potential issues of the new image product, which are of reference values for future studies. Second, this study uses the conceptual framework of recovery trajectory to quantify temporal patterns of NTL radiance in a disaster. This brings cutting-edge remote sensing techniques to resolve challenges in resilience assessment. The spatial analyses of disturbance (D NTL ) and recovery (R NTL ) revealed geographical disparities of human activities during the disaster, which provide novel insights into community recovery and resilience. We recognize the complexity of community resilience, which cannot be adequately measured using the two indicators. Instead, the study demonstrates the value of NTL images in providing timely and high-resolution data of human dynamics in disasters, which serves as an empirical grounding for developing and validating resilience measurement.

Interpretation of NTL Spatial Pattern
By comparing the disturbance (D NTL ) and recovery rate (R NTL ), the analysis revealed geographical disparities of recovery patterns in different CBSAs. In general, populated areas near the landfall location (including New York City) experienced a large disturbance but fast recovery, indicating the strong ability to bounce back in these areas. At the pixel level, the hurricane caused a large disturbance near the landfall location, including Ocean City in Maryland and Atlantic City in New Jersey. However, this region did not recover as fast as New York City. Several small CBSAs northeast to the landfall location (e.g., Springfield and Torrington in Connecticut and Pittsfield in Massachusetts) had a relatively smaller disturbance and recovery slower, which is possibly due to the impact of Nor'easter.
The analysis shows that the NTL metrics (D NTL and R NTL ) are significantly correlated (p < 0.05 and R 2 > 0.13) with 6 of the 36 variables, including wind speed, distance to hurricane, and ratios of keyword 'sandy', 'damage', and 'electric' in Twitter. Although the relatively low R 2 suggests limited prediction power of the regression models, the detected significant relations provide insights to the potential factors causing the NTL changes. The significant relations of the NTL metrics with wind speed indicate that the NTL fluctuation is mainly attributed to the hazardous weather condition brought by the hurricane. This finding confirmed that damaging wind is one of the primary causes of hurricane-inflicted loss of life and property damage as documented in the literature [34,[60][61][62]. Additionally, the high-wind areas are mostly located near the hurricane track, where the intensity of human activities were significantly reduced due to various reasons. For example, the largest NTL reduction is observed in the New York/New Jersey Bight area, which is a highly populated coastal area to the northeast of the landfall location. During the hurricane, a large population was evacuated from the predicted impact area, especially in New York City and surrounding cities in New Jersey [63]. Also, the hurricane has caused businesses and road closures, which lead to reduced visitations and traffic [64,65]. The evacuation and business closures are possibly the cause of NTL radiance reduction observed in urban areas. The relatively small NTL reduction in the center of New York City may reflect the local gathering of emergency responders [64]. The significant relation between NTL disturbance (D NTL ) and the ratio of keyword 'sandy' demonstrates a direct link between social media activity and NTL disturbance. Intensive discussions about the hurricane on Twitter occurred in areas that experienced large NTL disturbance. The literature shows that Twitter data is an effective indicator of human dynamics in disasters [54,66,67]. The correlation between NTL and Twitter further confirms the observed NTL changes likely reflect socio-economic impacts instead of other factors. Specifically, NTL disturbance (D NTL ) is significantly related to the ratio of keyword 'damage' in Twitter, which implies that the NTL radiance reduction may reflect damage to properties and infrastructures. Literature documented severe property damage in Ocean County in New Jersey [68], which is consistent with the large reduction of NTL observed in Figure 9a. The significant relation between the ratio of keyword 'electric' and NTL recovery (R NTL ) implies that areas concerned about electric issues recovered more quickly. This result agrees with the literature that restoration of electric supply played an important role in community and business recovery [69,70]. Additionally, D NTL has a significant relation with distance to hurricane track at the CBSA level, which confirms that the NTL disturbance is related to hazardous conditions in the hurricane. Areas closer to hurricane track were likely to experience strong wind, storm surge, and flooding brought by the hurricane, which lead to reduced intensity of human activity (e.g., evacuation, road/business closure).
Despite the detected significant relations, further research is needed to build robust models to explain the underlying causes of the observed NTL changes. The relatively low R 2 in the regression models may be attributed to insufficient high-quality images for continuous monitoring, as well as the lack of fine-resolution human dynamics data to compare with the NTL changes. Also, resampling the NTL data to the same scale of the exploratory variables may result in potential data loss or bias, thus affecting the model fitting. Future work is needed to ground-truthing the detected NTL patterns from other data sources, such as real-time traffic data [51], mobile phone tracking data [71], and survey data.
In general, larger disturbance (large absolute D NTL ) and faster recovery (large R NTL ) can be observed in urban land cover and roads. In contrast, most rural land cover and road types had smaller NTL changes (small absolute D NTL and small R NTL ) during the hurricane. This rural-urban disparity possibly indicates that NTL radiance is more sensitive to human dynamics in the populated area. Meanwhile, some rural areas in the state of Virginia, New York, and Connecticut still have not recovered to the pre-disaster condition on 17 November. The slow recovery potentially indicates low resilience in these rural areas. In addition to the large cities where public attention and resources are concentrated, more actions are needed to assist the recovery of rural areas that are often marginalized in disaster management [72,73].

Utility of the Black Marble Product in Disaster Monitoring
This study has attempted the use of new moonlight-adjusted NTL daily images to detect short-term recovery trajectories of human activities in Hurricane Sandy. Despite the limited high-quality images covering the study area during the hurricane, a V-shape pattern of NTL radiance has been detected in 82.1% of the studied CBSAs, which generally resemble the pattern of a recovery trajectory observed in other studies [33]. This result confirmed the utility of the daily Black Marble images in sensing disturbance and recovery of human activities in the hurricane. The remaining 12 CBSAs (17.9%) showing other patterns may be explained by a few possible reasons. First, these places were not impacted by the hurricane, and the NTL changes reflect normal activities (e.g., traffic congestion, public events, and power outage). Second, the larger in-disaster NTL can result from the evacuated population, which normally illuminated the evacuation sites and increase the local NTL. Third, inherent noise in Black Marble images, such as atmospheric reflection and lightning, may also introduce bias to the NTL radiance. Further investigations are needed to explain the real causes of the other patterns.
Despite the detected patterns of NTL changes, due to the characteristics of the Black Marble images, the analysis is still limited in the following aspects. First, the utility of daily NTL images is highly dependent on the atmospheric condition. Due to the cloudy sky during the hurricane, only images collected on some days have enough high-quality pixels covering the study area. Thus, the analysis is limited to three sampling points during the hurricane. Second, mosaicking images with different acquisition times may cause inconsistent NTL radiance in the study area. To overcome this issue, this study only selected days when the images were from a single scan, which ensures consistent NTL radiance over space. However, this selection criterion further reduces useable images for analysis. Ideally, the radiance inconsistency at different acquisition times can be calibrated before temporal analysis, which requires future work. Third, the influence of the other storm (Nor'easter Storm) affected the detected recovery pattern. Though dates around the Nor'easter are excluded in this study, NTL radiance at the post-disaster phase may be affected by the Nor'easter. The results interpreted in this study may reflect overlapped impacts and community responses in both disasters. In summary, these limitations indicate that the daily NTL images may not continuously provide high-quality data in hurricanes and flooding events when the atmospheric condition is not stable. Systematic data filtering, processing, and analysis are needed to develop valid analysis results. Instead, the daily NTL images may be more suited in monitoring disaster processes with a clear sky, such as an earthquake and tsunami. Although the moonlight-adjusted Black Marble Product brings advantages in resolution, data quality, and superior calibration algorithms, their applications need to be carefully designed according to specific disaster conditions.

Future Work
A few improvements in NTL products can be conducted in future studies. First, the resolution of current NTL data can be improved with additional data layers. Despite the 500 m (15 arc seconds) resolution of the Black Marble images is already an improvement from the older generation DMSP-OLS products, the pixel size is still too large to be precisely assigned to different land-use types or damaged structures, which creates difficulty to interpret the detected radiance changes. A possible solution is upsampling the NTL data to a finer resolution with additional data layers.  [25]. Second, new techniques can be applied to make use of low-quality pixels in NTL images. In this study, only high-quality pixels were used to analyze the daily NTL changes, leading to the limited coverage of the disaster area. In future studies, low-quality pixels may be calibrated by simulating atmospheric refraction and interpolating using neighboring high-quality pixels. Including low-quality pixels in analysis can significantly increase the monitoring area to obtain continuous data about human dynamics in disaster areas. Third, the validation issue needs to be solved. Currently, the detected NTL patterns cannot be fully explained by the selected variables, thus the causes of fluctuation remain uncertain. In the future, the NTL radiance will be compared with additional types of data (e.g., traffic data, mobile phone tracking data, and survey data) to fully understand the relation between NTL fluctuation and community recovery and resilience. Additionally, through applications in more disaster events, the utility of NTL data for resilience assessment can be further tested and the analytical method can be refined.

Conclusions
This study evaluated the use of NASA's new moonlight-adjusted Black Marble Products (SNPP-VIIRS VNP46A2) in detecting human activities during a major disaster (Hurricane Sandy). The study documented the analytical workflow from data collection, processing to data analysis and interpretation. The spatial analysis revealed the spatial variation of NTL radiance changes during the disaster. Univariate regression analyses were conducted to associate the NTL changes with auxiliary data, aiming to find explaining factors of the observed spatial variance. This study not only confirmed the merits of the NTL products in observing community recovery in disasters, but also identified issues that need to be resolved in future analysis. The findings and limitations discussed in this study recommend further research on associating additional ancillary datasets on environmental, socioeconomic, and demographic aspects. Overall, this study shows that daily NTL data is of great potential for exploring the impact on human activities during disasters and calls for more contributions in the methods and applications of daily SNPP-VIIRS NTL data. For broader implications, our framework on daily NTL data can be applied in other disasters and broader fields, such as detecting population movement, economic growth, and urbanization at both short-term and long-term scales.