Patterns in Land-use change and earthquake damage

Published:

Definitions

  • Land Use Change is defined as the number of observed changes in LANDSAT images (@NASADataLandsatScience) between 1 Jan 1988 and 1 Jan 2018, processed using the YATSM algorithm (@ZhuContinuouschangedetection2014)

  • Village lights are village-level averages of lights from the DMSP archive (@NOAANighttimeLightsAnnual)

Land Use Change and Population

Does the land-use change data correspond with observed population

#%%
change %<>% mutate(
    ln_hh = log(HHs2011),
    ln_popn = log(Popn2011),
    ln_popden = log(Popden2011),
    ln_chg    = log(chg_mean+1)
)
ggplot(data = change,aes(x=ln_popn, y = ln_chg)) +
    geom_point(aes(size=Popden2011, colour=ZONE_NAME),alpha=0.5) +
    geom_smooth()+labs(title='Log Population ~ Log Land-use Change')

plot of chunk unnamed-chunk-2

#%%
lm0 = robustify(felm(ln_chg ~ ln_popn, data=change))
stargazer(lm0,
  type='html'
  # type='text'
  )
Dependent variable:
ln_chg
ln_popn0.248***
(0.027)
Constant-0.913***
(0.228)
Observations627
R20.122
Adjusted R20.120
Residual Std. Error0.403 (df = 625)
Note:*p<0.1; **p<0.05; ***p<0.01
#%%

Yes.

Land use change and light

Does the land-use change data correspond with observed night-lights (which themselves are crude proxies of economic activity)?

#%%
change %>% left_join(lights, by=c('DIST_ID','VDC_ID')) -> light_luc
light_luc$ln_light = log(light_luc$village_light_mean+1)
lights$ln_light = log(lights$village_light_mean+1)
# ggplot() +
#     geom_density(data = light_luc,
#         mapping = aes(x=ln_light), colour='red')+
#     geom_density(data = lights,
#         mapping = aes(x=ln_light), colour='blue')
ggplot(data = light_luc,aes(x=ln_light, y = ln_chg)) +
    geom_point(aes(size=Popden2011, colour=ZONE_NAME.x),alpha=0.5) +
    geom_smooth(method='lm')+labs(title='Log Population ~ Log Light')

plot of chunk unnamed-chunk-5

#%%
lightreg = robustify(felm( ln_chg ~ ln_light, data=light_luc))
stargazer(lightreg,
  # type='text'
  type='html'
  )
Dependent variable:
ln_chg
ln_light0.291***
(0.020)
Constant1.085***
(0.018)
Observations627
R20.063
Adjusted R20.062
Residual Std. Error0.417 (df = 625)
Note:*p<0.1; **p<0.05; ***p<0.01
#%%

Yes.

Earthquake Damage

Calculate distance to epicenters of the two major earthquakes

Apr 25 - 27° 50′ 13.2″ N, 86° 4′ 37.2″ E

May 12 - 28° 13′ 48″ N, 84° 43′ 51.6″ E

# Earthquake epicenter
date = c('25 April 2015', '12 May 2015')
lat = c(28.23, 27.83)
lon = c(84.731, 86.077)
epc = data.frame(date, lat,lon)
epicenter = st_as_sf(epc, coords=c('lon','lat'), crs= 4326)

villages = st_read(paste0(data,'/processed/change_by_polygon.shp'), crs=4326)
## Reading layer `change_by_polygon' from data source `/home/alal/Desktop/code/0_research/eq-svy-eda/data/processed/change_by_polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 627 features and 16 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 84.42 ymin: 26.92 xmax: 86.69 ymax: 28.75
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
centroids = st_read(paste0(data,'/processed/village_centroids.shp'), crs=4326)
## Reading layer `village_centroids' from data source `/home/alal/Desktop/code/0_research/eq-svy-eda/data/processed/village_centroids.shp' using driver `ESRI Shapefile'
## Simple feature collection with 627 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 84.46 ymin: 26.97 xmax: 86.67 ymax: 28.64
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
#%%
tm_shape(villages) +
    tm_fill(col='chg_mean',
    palette = viridis(n=5, direction = -1,option = "C"))

plot of chunk unnamed-chunk-8

#%%

Centroids check

tm_shape(villages) +
    tm_borders() +
    tm_fill(col='Popn2011', palette = viridis(n=5, direction = 1,option = "C"))+
tm_shape(epicenter) +
    tm_dots(col="green",
            border.col="black",size=.3) +
tm_shape(centroids) +
    tm_dots(col='red')

plot of chunk unnamed-chunk-9

Calculate Haussdorf distance to epicenters

Pick the closer epicenter for each district

distances = as.data.frame(t(st_distance(epicenter, centroids)))
distances$mindist = apply(distances, 1, min)
centroids$disteq = distances$mindist

df = as.data.frame(centroids)
df %<>% select(-geometry)
villages %>% left_join(df, by=c('DIST_ID','VDC_ID')) -> vil_dists

Plot Haussdorf distance from Earthquake epicenter

plot(vil_dists['disteq'])

plot of chunk unnamed-chunk-11

#%%
#%%
dists = as.data.frame(vil_dists)
dists %<>% select(-geometry)
fwrite(dists,paste0(data,'/processed/villages_changes_distance.csv'))

Building damage

#%%
buildings = fread(paste0(data,'/interim/buildings.csv'))
## 
Read 15.2% of 1052948 rows
Read 34.2% of 1052948 rows
Read 57.0% of 1052948 rows
Read 80.7% of 1052948 rows
Read 1052948 rows and 125 (of 125) columns from 0.567 GB file in 00:00:07
#%% generate damage flags
buildings %<>% mutate(
    no_damage       = 1*(damage_overall_collapse == 'None' | damage_overall_collapse == ''),
    light_damage    = 1*(damage_overall_collapse == 'Insignificant/light'),
    moderate_damage = 1*(damage_overall_collapse == 'Moderate-Heavy'),
    severe_damage   = 1*(damage_overall_collapse == 'Severe-Extreme')
)

buildings %>% group_by(District,`VDC-Municipality-Code`) %>%
    select(no_damage, light_damage, moderate_damage, severe_damage, age_building) %>%
    summarize_all(funs(mean)) %>%
    rename(DIST_ID = District, VDC_ID=`VDC-Municipality-Code`) -> vlevel_buildings

vlevel_buildings %<>% inner_join(dists, by = c('DIST_ID', 'VDC_ID'))

Mapping damage

villages %>%
    left_join(vlevel_buildings, by=c('DIST_ID','VDC_ID')) ->
    vill_buildings
#%%
tm_shape(vill_buildings) +
    tm_borders() +
    tm_fill(col='severe_damage',palette = viridis(n=5, direction = -1,option = "C"))+
tm_shape(epicenter) +
    tm_dots(col="black",
        border.col="black",size=.3)

plot of chunk unnamed-chunk-15

#%%
#%%
tm_shape(vill_buildings) +
    tm_borders() +
    tm_fill(col='moderate_damage',palette = viridis(n=5, direction = -1,option = "C"))+
tm_shape(epicenter) +
    tm_dots(col="black",
        border.col="black",size=.3)

plot of chunk unnamed-chunk-16

#%%
vlevel_buildings %<>%
    mutate(lndist = log(disteq), ln_chg = log(chg_mean+1),
    ln_population = log(Popn2011)) %>% ungroup()

Conditional Expectation Functions

severe = binscatter('severe_damage ~ ln_chg + lndist + ln_population',
  key_var = 'ln_chg', data=vlevel_buildings, partial=T)
severe + labs(title='Severe earthquake damage ~ Land use change',
  subtitle='Partialling out distance to epicentre and population')

plot of chunk unnamed-chunk-18

moderate = binscatter('moderate_damage ~ ln_chg + lndist + ln_population',
  key_var = 'ln_chg', data=vlevel_buildings, partial=T)
moderate + labs(title='moderate earthquake damage ~ Land use change',
  subtitle='Partialling out distance to epicentre and population')

plot of chunk unnamed-chunk-18

light = binscatter('light_damage ~ ln_chg + lndist + ln_population',
  key_var = 'ln_chg', data=vlevel_buildings, partial=T)
light + labs(title='light earthquake damage ~ Land use change',
  subtitle='Partialling out distance to epicentre and population')

plot of chunk unnamed-chunk-18

Summary Statistics

#%%
svars <- vlevel_buildings %>% select(severe_damage,
    moderate_damage, light_damage, chg_mean, lndist) %>% ungroup()

knitr::kable(summary(svars))
 severe_damagemoderate_damagelight_damagechg_meanlndist
 Min. :0.0000Min. :0.0042Min. :0.0000Min. :0.00Min. : 7.73
 1st Qu.:0.06291st Qu.:0.11431st Qu.:0.02141st Qu.:1.871st Qu.:10.40
 Median :0.1291Median :0.2615Median :0.0717Median :2.36Median :10.79
 Mean :0.1591Mean :0.2705Mean :0.1248Mean :2.25Mean :10.66
 3rd Qu.:0.21723rd Qu.:0.40023rd Qu.:0.19663rd Qu.:2.813rd Qu.:11.07
 Max. :0.6569Max. :0.8153Max. :0.6620Max. :4.91Max. :11.54

Regressions

#%%
lm1 = robustify(felm(severe_damage   ~ ln_chg + ln_population + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))
lm2 = robustify(felm(moderate_damage ~ ln_chg + ln_population + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))
lm3 = robustify(felm(light_damage    ~ ln_chg + ln_population + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))

lm4 = robustify(felm(severe_damage   ~ ln_chg + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))
lm5 = robustify(felm(moderate_damage ~ ln_chg + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))
lm6 = robustify(felm(light_damage    ~ ln_chg + lndist|DIST_ID
    |0|DIST_ID, data=vlevel_buildings))
#%%
stargazer(lm4, lm1,lm5,lm2,lm6,lm3,
  covariate.labels= c('Log Land-use Change', 'Log Population',
  'Log Distance from Nearest Earthquake Epicenter'),
  add.lines = list(c("District Fixed Effects", "X","X","X","X","X","X")),
  # type='text'
  type='html'
  )
Dependent variable:
severe_damagemoderate_damagelight_damage
(1)(2)(3)(4)(5)(6)
Log Land-use Change0.0250.0410.044*0.074**0.005-0.011
(0.026)(0.029)(0.026)(0.030)(0.038)(0.034)
Log Population-0.021**-0.038*0.020**
(0.008)(0.021)(0.008)
Log Distance from Nearest Earthquake Epicenter0.0110.0160.103**0.111***0.050***0.046***
(0.019)(0.019)(0.044)(0.043)(0.015)(0.015)
District Fixed EffectsXXXXXX
Observations565565565565565565
R20.1230.1300.3700.3820.4700.475
Adjusted R20.1040.1100.3570.3670.4580.463
Residual Std. Error0.115 (df = 552)0.115 (df = 551)0.143 (df = 552)0.142 (df = 551)0.098 (df = 552)0.097 (df = 551)
Note:*p<0.1; **p<0.05; ***p<0.01

Upshot : Villages that urbanised more quickly suffered more damage in the earthquake after controlling from population and distance from the epicenters of the two major earthquakes and district level unobservables.

Bibliography