Spatial Analysis: 03
Tzai-Hung Wen
2025-03-10
Spatial Intersection: Point-In-Polygon
rm(list = ls())
library(sf)
library(tmap)
load("./data/Sample2.RData")1.1 Tornado Damage Assessment
Mapping
us_bg <- tm_shape(us_states_sf) + tm_polygons("grey90")
tornado <- tm_shape(torn_sf) + tm_dots(size = 0.04, shape = 1, fill_alpha = 0.5)
us_border <- tm_shape(us_states_sf) + tm_borders(col = "black") + tm_layout(frame = F)
us_bg + tornado + us_border
Checking CRS and data attribute
st_crs(us_states_sf)Coordinate Reference System:
User input: +proj=longlat +ellps=WGS84
wkt:
GEOGCRS["unknown",
DATUM["Unknown based on WGS84 ellipsoid",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]head(torn_sf)Simple feature collection with 6 features and 23 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -94.37 ymin: 30.67 xmax: -87.92 ymax: 34.48
Geodetic CRS: +proj=longlat +ellps=WGS84
TORNADX020 YEAR NUM STATE MONTH DAY DATE TOR_NO NO_STS STATE_TOR SEGNO STLAT STLON
0 1 1950 1 01 4 18 4/18/1950 53 1 1 1 30.67 88.20
1 2 1950 2 01 4 18 4/18/1950 54 1 1 1 30.70 87.92
2 3 1950 1 05 1 13 1/13/1950 4 1 1 1 34.40 94.37
3 4 1950 2 05 2 12 2/12/1950 19 1 1 1 34.48 92.40
4 5 1950 3 05 2 12 2/12/1950 23 1 1 1 33.27 92.95
5 6 1950 4 05 3 26 3/26/1950 33 1 1 1 34.12 93.07
SPLAT SPLON LGTH WIDTH FATAL INJ DAMAGE F_SCALE coords.x1 coords.x2 geometry
0 30.85 88.10 140 30 0 15 4 3 -88.20 30.67 POINT (-88.2 30.67)
1 0.00 0.00 20 45 0 0 3 2 -87.92 30.70 POINT (-87.92 30.7)
2 0.00 0.00 6 5 1 1 3 3 -94.37 34.40 POINT (-94.37 34.4)
3 0.00 0.00 1 30 0 0 3 2 -92.40 34.48 POINT (-92.4 34.48)
4 33.35 92.95 57 30 0 0 4 2 -92.95 33.27 POINT (-92.95 33.27)
5 34.32 92.88 174 45 0 3 4 2 -93.07 34.12 POINT (-93.07 34.12)summary(torn_sf) TORNADX020 YEAR NUM STATE MONTH
Min. : 1 Min. :1950 Min. : 1.00 48 : 7122 Min. : 1.00
1st Qu.:11734 1st Qu.:1970 1st Qu.: 7.00 40 : 2959 1st Qu.: 5.00
Median :23466 Median :1983 Median : 18.00 20 : 2905 Median : 6.00
Mean :23466 Mean :1982 Mean : 29.23 12 : 2788 Mean : 6.02
3rd Qu.:35199 3rd Qu.:1995 3rd Qu.: 38.00 31 : 2288 3rd Qu.: 7.00
Max. :46931 Max. :2004 Max. :232.00 19 : 1915 Max. :12.00
(Other):26954
DAY DATE TOR_NO NO_STS STATE_TOR
Min. : 1.0 4/3/1974 : 128 Min. : 1.0 Min. :1.000 Min. :0.0000
1st Qu.: 8.0 6/24/2003: 95 1st Qu.: 214.0 1st Qu.:1.000 1st Qu.:1.0000
Median :16.0 1/21/1999: 86 Median : 439.0 Median :1.000 Median :1.0000
Mean :15.7 5/30/2004: 81 Mean : 483.9 Mean :1.007 Mean :0.9932
3rd Qu.:23.0 5/4/2003 : 81 3rd Qu.: 695.0 3rd Qu.:1.000 3rd Qu.:1.0000
Max. :31.0 5/18/1995: 79 Max. :1817.0 Max. :3.000 Max. :1.0000
(Other) :46381
SEGNO STLAT STLON SPLAT SPLON
Min. :1 Min. :18.20 Min. : 64.90 Min. : 0.00 Min. : 0.00
1st Qu.:1 1st Qu.:33.13 1st Qu.: 86.85 1st Qu.: 0.00 1st Qu.: 0.00
Median :1 Median :37.08 Median : 94.23 Median : 0.00 Median : 0.00
Mean :1 Mean :37.14 Mean : 93.07 Mean :16.59 Mean : 40.96
3rd Qu.:1 3rd Qu.:41.07 3rd Qu.: 98.73 3rd Qu.:36.20 3rd Qu.: 91.38
Max. :1 Max. :61.02 Max. :163.53 Max. :61.02 Max. :163.53
LGTH WIDTH FATAL INJ
Min. : 0.00 Min. : 0.00 Min. : 0.00000 Min. : 0.000
1st Qu.: 1.00 1st Qu.: 3.00 1st Qu.: 0.00000 1st Qu.: 0.000
Median : 5.00 Median : 9.00 Median : 0.00000 Median : 0.000
Mean : 35.09 Mean : 23.55 Mean : 0.09987 Mean : 1.714
3rd Qu.: 27.00 3rd Qu.: 21.00 3rd Qu.: 0.00000 3rd Qu.: 0.000
Max. :3000.00 Max. :999.00 Max. :116.00000 Max. :1740.000
DAMAGE F_SCALE coords.x1 coords.x2
Min. : 0.000 Min. :-9.0000 Min. :-163.53 Min. :18.20
1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.: -98.73 1st Qu.:33.13
Median : 2.000 Median : 1.0000 Median : -94.23 Median :37.08
Mean : 2.271 Mean : 0.4819 Mean : -93.07 Mean :37.14
3rd Qu.: 4.000 3rd Qu.: 1.0000 3rd Qu.: -86.85 3rd Qu.:41.07
Max. :1000.000 Max. : 5.0000 Max. : -64.90 Max. :61.02
geometry
POINT :46931
epsg:NA : 0
+proj=long...: 0
st_geometry(torn_sf)Geometry set for 46931 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -163.53 ymin: 18.2 xmax: -64.9 ymax: 61.02
Geodetic CRS: +proj=longlat +ellps=WGS84
First 5 geometries:POINT (-88.2 30.67)
POINT (-87.92 30.7)
POINT (-94.37 34.4)
POINT (-92.4 34.48)
POINT (-92.95 33.27)Attribute selection
index <- us_states_sf$STATE_NAME == "Texas" | us_states_sf$STATE_NAME == "New Mexico" |
us_states_sf$STATE_NAME == "Oklahoma" | us_states_sf$STATE_NAME == "Arkansas"
AoI_sf <- us_states_sf[index,]
AoI_bg <- tm_shape(AoI_sf) + tm_polygons("grey90")
AoI_bg 
Clip
torn_clip_sf <- torn_sf[AoI_sf,]
torn_clip <- tm_shape(torn_clip_sf) + tm_dots(fill = "red", size = 0.04, fill_alpha = 0.5)
AoI_bg + torn_clip
head(torn_clip_sf)Simple feature collection with 6 features and 23 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -94.37 ymin: 33.27 xmax: -91.83 ymax: 36.15
Geodetic CRS: +proj=longlat +ellps=WGS84
TORNADX020 YEAR NUM STATE MONTH DAY DATE TOR_NO NO_STS STATE_TOR SEGNO STLAT STLON
2 3 1950 1 05 1 13 1/13/1950 4 1 1 1 34.40 94.37
3 4 1950 2 05 2 12 2/12/1950 19 1 1 1 34.48 92.40
4 5 1950 3 05 2 12 2/12/1950 23 1 1 1 33.27 92.95
5 6 1950 4 05 3 26 3/26/1950 33 1 1 1 34.12 93.07
6 7 1950 5 05 3 26 3/26/1950 34 1 1 1 36.15 91.83
7 8 1950 6 05 3 26 3/26/1950 35 1 1 1 34.70 92.35
SPLAT SPLON LGTH WIDTH FATAL INJ DAMAGE F_SCALE coords.x1 coords.x2 geometry
2 0.00 0.00 6 5 1 1 3 3 -94.37 34.40 POINT (-94.37 34.4)
3 0.00 0.00 1 30 0 0 3 2 -92.40 34.48 POINT (-92.4 34.48)
4 33.35 92.95 57 30 0 0 4 2 -92.95 33.27 POINT (-92.95 33.27)
5 34.32 92.88 174 45 0 3 4 2 -93.07 34.12 POINT (-93.07 34.12)
6 36.20 91.75 57 60 0 1 5 3 -91.83 36.15 POINT (-91.83 36.15)
7 34.80 92.22 104 180 0 7 5 2 -92.35 34.70 POINT (-92.35 34.7)Point-In-Polygon: Using st_join()
head(AoI_torn_sf)
Simple feature collection with 6 features and 74 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -94.37 ymin: 33.27 xmax: -91.83 ymax: 36.15
Geodetic CRS: +proj=longlat +ellps=WGS84
TORNADX020 YEAR NUM STATE MONTH DAY DATE TOR_NO NO_STS STATE_TOR SEGNO STLAT STLON
2 3 1950 1 05 1 13 1/13/1950 4 1 1 1 34.40 94.37
3 4 1950 2 05 2 12 2/12/1950 19 1 1 1 34.48 92.40
4 5 1950 3 05 2 12 2/12/1950 23 1 1 1 33.27 92.95
5 6 1950 4 05 3 26 3/26/1950 33 1 1 1 34.12 93.07
6 7 1950 5 05 3 26 3/26/1950 34 1 1 1 36.15 91.83
7 8 1950 6 05 3 26 3/26/1950 35 1 1 1 34.70 92.35
SPLAT SPLON LGTH WIDTH FATAL INJ DAMAGE F_SCALE coords.x1 coords.x2 AREA STATE_NAME
2 0.00 0.00 6 5 1 1 3 3 -94.37 34.40 52912.79 Arkansas
3 0.00 0.00 1 30 0 0 3 2 -92.40 34.48 52912.79 Arkansas
4 33.35 92.95 57 30 0 0 4 2 -92.95 33.27 52912.79 Arkansas
5 34.32 92.88 174 45 0 3 4 2 -93.07 34.12 52912.79 Arkansas
6 36.20 91.75 57 60 0 1 5 3 -91.83 36.15 52912.79 Arkansas
7 34.80 92.22 104 180 0 7 5 2 -92.35 34.70 52912.79 Arkansas
STATE_FIPS SUB_REGION STATE_ABBR POP1990 POP1997 POP90_SQMI HOUSEHOLDS MALES FEMALES
2 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
3 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
4 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
5 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
6 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
7 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
WHITE BLACK AMERI_ES ASIAN_PI OTHER HISPANIC AGE_UNDER5 AGE_5_17 AGE_18_29 AGE_30_49
2 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
3 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
4 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
5 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
6 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
7 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
AGE_50_64 AGE_65_UP NEVERMARRY MARRIED SEPARATED WIDOWED DIVORCED HSEHLD_1_M HSEHLD_1_F
2 324647 350058 379361 1094739 35256 164342 161212 80266 133512
3 324647 350058 379361 1094739 35256 164342 161212 80266 133512
4 324647 350058 379361 1094739 35256 164342 161212 80266 133512
5 324647 350058 379361 1094739 35256 164342 161212 80266 133512
6 324647 350058 379361 1094739 35256 164342 161212 80266 133512
7 324647 350058 379361 1094739 35256 164342 161212 80266 133512
MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD HSE_UNITS VACANT OWNER_OCC RENTER_OCC MEDIAN_VAL
2 248131 279227 13819 68471 1000667 109488 619938 271241 46300
3 248131 279227 13819 68471 1000667 109488 619938 271241 46300
4 248131 279227 13819 68471 1000667 109488 619938 271241 46300
5 248131 279227 13819 68471 1000667 109488 619938 271241 46300
6 248131 279227 13819 68471 1000667 109488 619938 271241 46300
7 248131 279227 13819 68471 1000667 109488 619938 271241 46300
MEDIANRENT UNITS_1DET UNITS_1ATT UNITS2 UNITS3_9 UNITS10_49 UNITS50_UP MOBILEHOME
2 230 708751 18175 31930 55914 35714 8740 131542
3 230 708751 18175 31930 55914 35714 8740 131542
4 230 708751 18175 31930 55914 35714 8740 131542
5 230 708751 18175 31930 55914 35714 8740 131542
6 230 708751 18175 31930 55914 35714 8740 131542
7 230 708751 18175 31930 55914 35714 8740 131542
NO_FARMS87 AVG_SIZE87 CROP_ACR87 AVG_SALE87 geometry
2 48242 298 9950401 68825 POINT (-94.37 34.4)
3 48242 298 9950401 68825 POINT (-92.4 34.48)
4 48242 298 9950401 68825 POINT (-92.95 33.27)
5 48242 298 9950401 68825 POINT (-93.07 34.12)
6 48242 298 9950401 68825 POINT (-91.83 36.15)
7 48242 298 9950401 68825 POINT (-92.35 34.7)Creating new data frame
AoI_torn_df$STATE_NAME <- droplevels(AoI_torn_df$STATE_NAME)
newdf<- data.frame(Name = AoI_torn_df$STATE_NAME,Damage= AoI_torn_sf$DAMAGE)
head(newdf)newdf$Damage <- as.integer(newdf$Damage)Crosstab analysis
count<- table(newdf$Damage, newdf$Name)
count
Arkansas New Mexico Oklahoma Texas
0 563 283 1280 4053
1 23 35 211 378
2 39 33 124 235
3 132 49 376 676
4 235 49 506 919
5 222 25 359 560
6 60 2 92 165
7 15 1 17 36
8 1 0 3 7
9 0 0 0 1
10 0 0 2 1
12 1 0 0 1
13 0 0 1 0
15 0 0 0 1
19 0 0 0 1
20 0 0 0 1
40 0 0 0 1
60 0 0 1 0
70 0 0 0 1
75 0 0 0 1
100 0 0 1 0
125 0 0 0 1
150 0 0 1 0
305 0 0 0 1
370 0 0 1 0
1000 0 0 1 0count2<- xtabs(~ newdf$Damage + newdf$Nam )
count2 newdf$Nam
newdf$Damage Arkansas New Mexico Oklahoma Texas
0 563 283 1280 4053
1 23 35 211 378
2 39 33 124 235
3 132 49 376 676
4 235 49 506 919
5 222 25 359 560
6 60 2 92 165
7 15 1 17 36
8 1 0 3 7
9 0 0 0 1
10 0 0 2 1
12 1 0 0 1
13 0 0 1 0
15 0 0 0 1
19 0 0 0 1
20 0 0 0 1
40 0 0 0 1
60 0 0 1 0
70 0 0 0 1
75 0 0 0 1
100 0 0 1 0
125 0 0 0 1
150 0 0 1 0
305 0 0 0 1
370 0 0 1 0
1000 0 0 1 01.2 Mapping Tornado Density
# (Counting Points in Polygons)
points_sf_joined <- st_join(torn_sf, us_states_sf, join = st_within)
head(points_sf_joined)Simple feature collection with 6 features and 77 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -94.37 ymin: 30.67 xmax: -87.92 ymax: 34.48
Geodetic CRS: +proj=longlat +ellps=WGS84
TORNADX020 YEAR NUM STATE MONTH DAY DATE TOR_NO NO_STS STATE_TOR SEGNO STLAT STLON
0 1 1950 1 01 4 18 4/18/1950 53 1 1 1 30.67 88.20
1 2 1950 2 01 4 18 4/18/1950 54 1 1 1 30.70 87.92
2 3 1950 1 05 1 13 1/13/1950 4 1 1 1 34.40 94.37
3 4 1950 2 05 2 12 2/12/1950 19 1 1 1 34.48 92.40
4 5 1950 3 05 2 12 2/12/1950 23 1 1 1 33.27 92.95
5 6 1950 4 05 3 26 3/26/1950 33 1 1 1 34.12 93.07
SPLAT SPLON LGTH WIDTH FATAL INJ DAMAGE F_SCALE coords.x1 coords.x2 AREA STATE_NAME
0 30.85 88.10 140 30 0 15 4 3 -88.20 30.67 51715.66 Alabama
1 0.00 0.00 20 45 0 0 3 2 -87.92 30.70 51715.66 Alabama
2 0.00 0.00 6 5 1 1 3 3 -94.37 34.40 52912.79 Arkansas
3 0.00 0.00 1 30 0 0 3 2 -92.40 34.48 52912.79 Arkansas
4 33.35 92.95 57 30 0 0 4 2 -92.95 33.27 52912.79 Arkansas
5 34.32 92.88 174 45 0 3 4 2 -93.07 34.12 52912.79 Arkansas
STATE_FIPS SUB_REGION STATE_ABBR POP1990 POP1997 POP90_SQMI HOUSEHOLDS MALES FEMALES
0 01 E S Cen AL 4040587 4298715 78 1506790 1936162 2104425
1 01 E S Cen AL 4040587 4298715 78 1506790 1936162 2104425
2 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
3 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
4 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
5 05 W S Cen AR 2350725 2529864 44 891179 1133076 1217649
WHITE BLACK AMERI_ES ASIAN_PI OTHER HISPANIC AGE_UNDER5 AGE_5_17 AGE_18_29 AGE_30_49
0 2975797 1020705 16506 21797 5782 24629 283295 775493 762897 1137367
1 2975797 1020705 16506 21797 5782 24629 283295 775493 762897 1137367
2 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
3 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
4 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
5 1944744 373912 12773 12530 6766 19876 164667 456464 416092 638797
AGE_50_64 AGE_65_UP NEVERMARRY MARRIED SEPARATED WIDOWED DIVORCED HSEHLD_1_M HSEHLD_1_F
0 558546 522989 754868 1791644 68002 276267 273511 138220 219858
1 558546 522989 754868 1791644 68002 276267 273511 138220 219858
2 324647 350058 379361 1094739 35256 164342 161212 80266 133512
3 324647 350058 379361 1094739 35256 164342 161212 80266 133512
4 324647 350058 379361 1094739 35256 164342 161212 80266 133512
5 324647 350058 379361 1094739 35256 164342 161212 80266 133512
MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD HSE_UNITS VACANT OWNER_OCC RENTER_OCC MEDIAN_VAL
0 417950 440377 21736 132896 1670379 163589 1061897 444893 53700
1 417950 440377 21736 132896 1670379 163589 1061897 444893 53700
2 248131 279227 13819 68471 1000667 109488 619938 271241 46300
3 248131 279227 13819 68471 1000667 109488 619938 271241 46300
4 248131 279227 13819 68471 1000667 109488 619938 271241 46300
5 248131 279227 13819 68471 1000667 109488 619938 271241 46300
MEDIANRENT UNITS_1DET UNITS_1ATT UNITS2 UNITS3_9 UNITS10_49 UNITS50_UP MOBILEHOME
0 229 1133927 31943 42295 120222 81983 20479 224307
1 229 1133927 31943 42295 120222 81983 20479 224307
2 230 708751 18175 31930 55914 35714 8740 131542
3 230 708751 18175 31930 55914 35714 8740 131542
4 230 708751 18175 31930 55914 35714 8740 131542
5 230 708751 18175 31930 55914 35714 8740 131542
NO_FARMS87 AVG_SIZE87 CROP_ACR87 AVG_SALE87 Counts AREA1 Density
0 43318 211 4496607 44053 1266 134018.9 [km^2] 0.009446430 [1/km^2]
1 43318 211 4496607 44053 1266 134018.9 [km^2] 0.009446430 [1/km^2]
2 48242 298 9950401 68825 1291 137058.1 [km^2] 0.009419366 [1/km^2]
3 48242 298 9950401 68825 1291 137058.1 [km^2] 0.009419366 [1/km^2]
4 48242 298 9950401 68825 1291 137058.1 [km^2] 0.009419366 [1/km^2]
5 48242 298 9950401 68825 1291 137058.1 [km^2] 0.009419366 [1/km^2]
geometry
0 POINT (-88.2 30.67)
1 POINT (-87.92 30.7)
2 POINT (-94.37 34.4)
3 POINT (-92.4 34.48)
4 POINT (-92.95 33.27)
5 POINT (-93.07 34.12)table1 <- table(points_sf_joined$STATE_NAME)
count_sf <- as.data.frame(table1)
colnames(count_sf) <- c("STATE_NAME","Counts")
head(count_sf)
library(dplyr)
us_states_sf<- left_join(us_states_sf, count_sf)Joining with `by = join_by(STATE_NAME, Counts)`head(us_states_sf)Simple feature collection with 6 features and 54 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -124.7314 ymin: 40.9943 xmax: -66.96985 ymax: 49
Geodetic CRS: +proj=longlat +ellps=WGS84
AREA STATE_NAME STATE_FIPS SUB_REGION STATE_ABBR POP1990 POP1997 POP90_SQMI
1 67286.88 Washington 53 Pacific WA 4866692 5604260 72
2 147236.03 Montana 30 Mtn MT 799065 888723 5
3 32161.66 Maine 23 N Eng ME 1227928 1244828 38
4 70810.15 North Dakota 38 W N Cen ND 638800 644782 9
5 77193.62 South Dakota 46 W N Cen SD 696004 736549 9
6 97799.49 Wyoming 56 Mtn WY 453588 484529 5
HOUSEHOLDS MALES FEMALES WHITE BLACK AMERI_ES ASIAN_PI OTHER HISPANIC AGE_UNDER5
1 1872431 2413747 2452945 4308937 149801 81483 210958 115513 214570 366780
2 306163 395769 403296 741111 2381 47679 4259 3635 12174 59257
3 465312 597850 630078 1208360 5138 5998 6683 1749 6829 85722
4 240878 318201 320599 604142 3524 25917 3462 1755 4665 47845
5 259034 342498 353506 637515 3258 50575 3123 1533 5252 54504
6 168839 227007 226581 427061 3606 9479 2806 10636 25751 34780
AGE_5_17 AGE_18_29 AGE_30_49 AGE_50_64 AGE_65_UP NEVERMARRY MARRIED SEPARATED WIDOWED
1 894607 900361 1531803 597853 575288 942004 2147036 71269 229431
2 162847 126122 238848 105494 106497 136569 365708 8132 44216
3 223280 222545 368850 164158 163373 233006 561930 12873 73288
4 127540 118007 174516 79837 91055 126912 292761 3769 37188
5 143958 121256 185013 88942 102331 128748 313946 5157 43070
6 100745 74654 140493 55721 47195 73555 207819 3809 20546
DIVORCED HSEHLD_1_M HSEHLD_1_F MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD HSE_UNITS VACANT
1 401417 212023 264297 490263 539004 35545 125022 2032378 159947
2 56907 36000 44491 85633 90893 5371 19416 361155 54992
3 88024 42075 66399 129896 140669 7884 29843 587045 121733
4 29473 27313 36640 71762 70612 3055 12217 276340 35462
5 36347 27680 40628 75798 76721 3927 15191 292436 33402
6 33545 19671 21616 53125 47675 3225 10710 203411 34572
OWNER_OCC RENTER_OCC MEDIAN_VAL MEDIANRENT UNITS_1DET UNITS_1ATT UNITS2 UNITS3_9
1 1171580 700851 93400 383 1272721 48086 61626 168162
2 205899 100264 56600 251 237533 8432 13640 26063
3 327888 137424 87400 358 378413 11753 36416 66214
4 157950 82928 50800 266 172938 10286 8748 24390
5 171161 87873 45200 242 202166 5249 8758 20411
6 114544 54295 61600 270 129197 6212 5911 15818
UNITS10_49 UNITS50_UP MOBILEHOME NO_FARMS87 AVG_SIZE87 CROP_ACR87 AVG_SALE87 Counts
1 218722 55864 187533 33559 480 8168454 87000 80
2 14174 2757 54021 24568 2451 17829766 62980 343
3 20813 5417 54532 6269 214 592309 64681 87
4 27427 2935 27055 35289 1143 28208099 62007 1128
5 19751 1891 31357 36376 1214 19641972 74761 1445
6 9501 917 33474 9205 3650 2838627 73517 548
AREA1 Density geometry
1 173800.40 [km^2] 0.0004602981 [1/km^2] MULTIPOLYGON (((-122.4007 4...
2 380316.56 [km^2] 0.0009018803 [1/km^2] MULTIPOLYGON (((-111.4746 4...
3 83108.37 [km^2] 0.0010468259 [1/km^2] MULTIPOLYGON (((-69.77779 4...
4 182881.38 [km^2] 0.0061679325 [1/km^2] MULTIPOLYGON (((-98.73006 4...
5 199514.44 [km^2] 0.0072425837 [1/km^2] MULTIPOLYGON (((-102.7879 4...
6 252857.54 [km^2] 0.0021672282 [1/km^2] MULTIPOLYGON (((-104.0531 4...area1<-st_area(us_states_sf) # unit: foot
library(units)
area1<-set_units(area1, km^2)
us_states_sf$AREA1 <- area1
us_states_sf$Density <- us_states_sf$Counts / us_states_sf$AREA1
head(us_states_sf)Simple feature collection with 6 features and 54 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -124.7314 ymin: 40.9943 xmax: -66.96985 ymax: 49
Geodetic CRS: +proj=longlat +ellps=WGS84
AREA STATE_NAME STATE_FIPS SUB_REGION STATE_ABBR POP1990 POP1997 POP90_SQMI
1 67286.88 Washington 53 Pacific WA 4866692 5604260 72
2 147236.03 Montana 30 Mtn MT 799065 888723 5
3 32161.66 Maine 23 N Eng ME 1227928 1244828 38
4 70810.15 North Dakota 38 W N Cen ND 638800 644782 9
5 77193.62 South Dakota 46 W N Cen SD 696004 736549 9
6 97799.49 Wyoming 56 Mtn WY 453588 484529 5
HOUSEHOLDS MALES FEMALES WHITE BLACK AMERI_ES ASIAN_PI OTHER HISPANIC AGE_UNDER5
1 1872431 2413747 2452945 4308937 149801 81483 210958 115513 214570 366780
2 306163 395769 403296 741111 2381 47679 4259 3635 12174 59257
3 465312 597850 630078 1208360 5138 5998 6683 1749 6829 85722
4 240878 318201 320599 604142 3524 25917 3462 1755 4665 47845
5 259034 342498 353506 637515 3258 50575 3123 1533 5252 54504
6 168839 227007 226581 427061 3606 9479 2806 10636 25751 34780
AGE_5_17 AGE_18_29 AGE_30_49 AGE_50_64 AGE_65_UP NEVERMARRY MARRIED SEPARATED WIDOWED
1 894607 900361 1531803 597853 575288 942004 2147036 71269 229431
2 162847 126122 238848 105494 106497 136569 365708 8132 44216
3 223280 222545 368850 164158 163373 233006 561930 12873 73288
4 127540 118007 174516 79837 91055 126912 292761 3769 37188
5 143958 121256 185013 88942 102331 128748 313946 5157 43070
6 100745 74654 140493 55721 47195 73555 207819 3809 20546
DIVORCED HSEHLD_1_M HSEHLD_1_F MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD HSE_UNITS VACANT
1 401417 212023 264297 490263 539004 35545 125022 2032378 159947
2 56907 36000 44491 85633 90893 5371 19416 361155 54992
3 88024 42075 66399 129896 140669 7884 29843 587045 121733
4 29473 27313 36640 71762 70612 3055 12217 276340 35462
5 36347 27680 40628 75798 76721 3927 15191 292436 33402
6 33545 19671 21616 53125 47675 3225 10710 203411 34572
OWNER_OCC RENTER_OCC MEDIAN_VAL MEDIANRENT UNITS_1DET UNITS_1ATT UNITS2 UNITS3_9
1 1171580 700851 93400 383 1272721 48086 61626 168162
2 205899 100264 56600 251 237533 8432 13640 26063
3 327888 137424 87400 358 378413 11753 36416 66214
4 157950 82928 50800 266 172938 10286 8748 24390
5 171161 87873 45200 242 202166 5249 8758 20411
6 114544 54295 61600 270 129197 6212 5911 15818
UNITS10_49 UNITS50_UP MOBILEHOME NO_FARMS87 AVG_SIZE87 CROP_ACR87 AVG_SALE87 Counts
1 218722 55864 187533 33559 480 8168454 87000 80
2 14174 2757 54021 24568 2451 17829766 62980 343
3 20813 5417 54532 6269 214 592309 64681 87
4 27427 2935 27055 35289 1143 28208099 62007 1128
5 19751 1891 31357 36376 1214 19641972 74761 1445
6 9501 917 33474 9205 3650 2838627 73517 548
AREA1 Density geometry
1 173800.40 [km^2] 0.0004602981 [1/km^2] MULTIPOLYGON (((-122.4007 4...
2 380316.56 [km^2] 0.0009018803 [1/km^2] MULTIPOLYGON (((-111.4746 4...
3 83108.37 [km^2] 0.0010468259 [1/km^2] MULTIPOLYGON (((-69.77779 4...
4 182881.38 [km^2] 0.0061679325 [1/km^2] MULTIPOLYGON (((-98.73006 4...
5 199514.44 [km^2] 0.0072425837 [1/km^2] MULTIPOLYGON (((-102.7879 4...
6 252857.54 [km^2] 0.0021672282 [1/km^2] MULTIPOLYGON (((-104.0531 4...library(pals) # collection of color palettes
plot(us_states_sf["Density"], breaks = "jenks", nbreaks = 6, pal=brewer.reds(6))
Buffer zone
st_crs(places_sf)Coordinate Reference System:
User input: +proj=lcc +datum=NAD27 +lon_0=-72d45 +lat_1=41d52 +lat_2=41d12 +lat_0=40d50 +x_0=182880.3657607315 +y_0=0 +units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["North American Datum 1927",
ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
LENGTHUNIT["metre",1]],
ID["EPSG",6267]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["unknown",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",40.8333333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",-72.75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",41.8666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",41.2,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",600000,
LENGTHUNIT["US survey foot",0.304800609601219],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",0,
LENGTHUNIT["US survey foot",0.304800609601219],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["US survey foot",0.304800609601219,
ID["EPSG",9003]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["US survey foot",0.304800609601219,
ID["EPSG",9003]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["unknown to WGS84",
METHOD["NTv2",
ID["EPSG",9615]],
PARAMETERFILE["Latitude and longitude difference file","@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat",
ID["EPSG",8656]]]]st_crs(places_sf)<- st_crs(blocks_sf)
head(places_sf)Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 539618.6 ymin: 164850.7 xmax: 561995.3 ymax: 181792.2
Projected CRS: +proj=lcc +datum=NAD27 +lon_0=-72d45 +lat_1=41d52 +lat_2=41d12 +lat_0=40d50 +x_0=182880.3657607315 +y_0=0 +units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
label geometry
1 NEW HAVEN POINT (549485.9 171645.5)
2 Westville POINT (539618.6 181792.2)
3 Fair Haven POINT (559858 174114.6)
4 New Haven POINT (550926.8 173018.1)
5 Fair Haven East POINT (561995.3 169399.7)
6 City Point POINT (547985.7 164850.7)block_bg <- tm_shape(blocks_sf) + tm_polygons("grey90")
places_pts <- tm_shape(places_sf) + tm_dots(fill = "red", size = 0.5)
block_bg + places_pts
places_buf_sf <- st_buffer(places_sf, dist = 3000)
summary(places_buf_sf) label geometry
Length:9 POLYGON :9
Class :character epsg:NA :0
Mode :character +proj=lcc ...:0 places_buf <- tm_shape(places_buf_sf) + tm_polygons("yellow")
block_bg + places_buf + places_pts
places_bufU_sf <- st_union(places_buf_sf)
head(places_bufU_sf)Geometry set for 1 feature
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 536618.6 ymin: 146437.8 xmax: 564995.3 ymax: 184792.2
Projected CRS: +proj=lcc +datum=NAD27 +lon_0=-72d45 +lat_1=41d52 +lat_2=41d12 +lat_0=40d50 +x_0=182880.3657607315 +y_0=0 +units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.datMULTIPOLYGON (((560838.8 148968.5, 560810.2 148...places_bufU <- tm_shape(places_bufU_sf) + tm_polygons("yellow")
block_bg + places_bufU + places_pts
Distance Analysis
checking attribute table
head(blocks_sf)Simple feature collection with 6 features and 28 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 534687.6 ymin: 177306.3 xmax: 569625.3 ymax: 188464.6
Projected CRS: +proj=lcc +datum=NAD27 +lon_0=-72d45 +lat_1=41d52 +lat_2=41d12 +lat_0=40d50 +x_0=182880.3657607315 +y_0=0 +units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
NEWH075H_ NEWH075H_I HSE_UNITS OCCUPIED VACANT P_VACANT P_OWNEROCC P_RENTROCC NEWH075P_
0 2 69 763 725 38 4.980341 0.393185 94.626474 2
1 3 72 510 480 30 5.882353 20.392157 73.725490 3
2 4 64 389 362 27 6.940874 57.840617 35.218509 4
3 5 68 429 397 32 7.459207 19.813520 72.727273 5
4 6 67 443 385 58 13.092551 80.361174 6.546275 6
5 7 133 588 548 40 6.802721 52.551020 40.646259 7
NEWH075P_I POP1990 P_MALES P_FEMALES P_WHITE P_BLACK P_AMERI_ES P_ASIAN_PI P_OTHER
0 380 2396 40.02504 59.97496 7.095159 87.020033 0.584307 0.041736 5.258765
1 385 3071 39.07522 60.92478 87.105177 10.452621 0.195376 0.521003 1.725822
2 394 996 47.38956 52.61044 32.931727 66.265060 0.100402 0.200803 0.502008
3 399 1336 42.66467 57.33533 11.452096 85.553892 0.523952 0.523952 1.946108
4 404 915 46.22951 53.77049 73.442623 24.371585 0.327869 1.420765 0.437158
5 406 1318 50.91047 49.08953 87.784522 7.435508 0.758725 0.834598 3.186646
P_UNDER5 P_5_13 P_14_17 P_18_24 P_25_34 P_35_44 P_45_54 P_55_64 P_65_74
0 12.813022 24.707846 7.888147 12.479132 16.026711 8.555927 5.759599 4.924875 4.048414
1 1.921198 2.474764 0.814067 71.149463 7.359166 4.037773 1.595571 1.758385 3.712146
2 10.441767 13.554217 5.722892 8.835341 17.670683 17.871486 8.734940 5.923695 7.931727
3 10.853293 17.739521 7.709581 12.425150 18.113772 10.853293 9.056886 6.287425 4.266467
4 6.229508 8.633880 2.950820 7.103825 17.267760 16.830601 8.415301 7.431694 14.426230
5 8.725341 8.194234 3.641882 10.091047 29.286798 12.898331 7.814871 7.814871 6.904401
P_75_UP geometry
0 2.796327 POLYGON ((540989.5 186028.3...
1 5.177467 POLYGON ((539949.9 187487.6...
2 3.313253 POLYGON ((537497.6 184616.7...
3 2.694611 POLYGON ((537497.6 184616.7...
4 10.710383 POLYGON ((536589.3 184217.5...
5 4.628225 POLYGON ((568032.4 183170.2...using st_centroid()
blocks_center_sf <- st_centroid(blocks_sf)è¦å‘Šï¼š st_centroid assumes attributes are constant over geometriesblocks_center <- tm_shape(blocks_center_sf) + tm_dots(fill = "blue", size = 0.2)
block_bg +blocks_center
blocks_center + places_pts
using st_distance()
distance.matrix <- st_distance(blocks_center_sf, places_sf)
library(units)
distance.matrix<-set_units(distance.matrix, km)using apply()
near_dist<- apply(distance.matrix, 1, mean)
near_dist[1][1] 6.979727xid<- which.min(near_dist)
near_dist[xid][1] 3.321075blocks_sf[xid,]Simple feature collection with 1 feature and 28 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 550067.1 ymin: 163375.2 xmax: 557408.4 ymax: 170921
Projected CRS: +proj=lcc +datum=NAD27 +lon_0=-72d45 +lat_1=41d52 +lat_2=41d12 +lat_0=40d50 +x_0=182880.3657607315 +y_0=0 +units=us-ft +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
NEWH075H_ NEWH075H_I HSE_UNITS OCCUPIED VACANT P_VACANT P_OWNEROCC P_RENTROCC NEWH075P_
109 111 26 807 774 33 4.089219 2.478315 93.43247 111
NEWH075P_I POP1990 P_MALES P_FEMALES P_WHITE P_BLACK P_AMERI_ES P_ASIAN_PI P_OTHER
109 545 1676 41.05012 58.94988 26.07399 63.90215 0.238663 0.059666 9.725537
P_UNDER5 P_5_13 P_14_17 P_18_24 P_25_34 P_35_44 P_45_54 P_55_64 P_65_74 P_75_UP
109 8.412888 15.99045 5.48926 10.02387 14.6778 8.353222 5.369928 5.548926 8.770883 17.36277
geometry
109 POLYGON ((550592.3 169501, ...blocks_sf$POP1990[xid][1] 1676min_dist<- apply(distance.matrix, 1, min)
sel_blocks<- min_dist < 1 # < 1 km
sel_sf <- blocks_sf[sel_blocks,]
sel_map <- tm_shape(sel_sf) + tm_polygons("yellow")
block_bg + sel_map + places_pts
using st_is_within_distance()
d1<- set_units(1, km)
sel_blocks = st_is_within_distance(blocks_center_sf, places_sf, dist = d1)
class(sel_blocks)[1] "sgbp" "list"sel_blocks<- lengths(sel_blocks) > 0
sel_sf <- blocks_sf[sel_blocks,]
sel_map <- tm_shape(sel_sf) + tm_polygons("yellow")
block_bg + sel_map + places_pts
choose specific one place
sel_blocks = st_is_within_distance(blocks_center_sf, places_sf[3,], dist = d1)
sel_blocks <- lengths(sel_blocks) > 0
sel_sf <- blocks_sf[sel_blocks,]
sel_map <- tm_shape(sel_sf) + tm_polygons("yellow")
block_bg + sel_map + places_pts
TOTAL_POP <- sum(sel_sf$POP1990)
TOTAL_WHITE <- sum(sel_sf$POP1990 * (sel_sf$P_WHITE/100) )accessibility analysis: using apply()
eth<- as.data.frame(blocks_center_sf[,14:18])
eth<- as.matrix(eth[,1:5])
eth<- apply(eth, 2, function(x) (x*blocks_center_sf$POP1990))
colnames(eth)<- c("White","Black","Native","Asian","Other")
eth.access <- xtabs(eth~sel_blocks)
eth.access <- as.data.frame(eth.access)
colnames(eth.access) <- c("Access","Eth5","Pop")
library(ggplot2)
ggplot(eth.access) + aes(fill=Access, y=Pop, x=Eth5) +
geom_bar(position="stack", stat="identity")
LS0tDQp0aXRsZTogIlNwYXRpYWwgQW5hbHlzaXM6IDAzIg0KYXV0aG9yOiAiVHphaS1IdW5nIFdlbiINCmRhdGU6ICcyMDI1LTAzLTEwJw0Kb3V0cHV0Og0KICBodG1sX25vdGVib29rOg0KICAgIHRvYzogdHJ1ZQ0KICAgIHRvY19kZXB0aDogNg0KICAgIHRvY19mbG9hdDogdHJ1ZQ0KLS0tDQoNCiMgU3BhdGlhbCBJbnRlcnNlY3Rpb246IFBvaW50LUluLVBvbHlnb24NCg0KYGBge3J9DQpybShsaXN0ID0gbHMoKSkNCg0KbGlicmFyeShzZikNCmxpYnJhcnkodG1hcCkNCg0KbG9hZCgiLi9kYXRhL1NhbXBsZTIuUkRhdGEiKQ0KYGBgDQoNCiMjIDEuMSBUb3JuYWRvIERhbWFnZSBBc3Nlc3NtZW50IA0KDQojIyMgTWFwcGluZw0KYGBge3J9DQp1c19iZyA8LSB0bV9zaGFwZSh1c19zdGF0ZXNfc2YpICsgdG1fcG9seWdvbnMoImdyZXk5MCIpDQp0b3JuYWRvIDwtIHRtX3NoYXBlKHRvcm5fc2YpICsgdG1fZG90cyhzaXplID0gMC4wNCwgc2hhcGUgPSAxLCBmaWxsX2FscGhhID0gMC41KQ0KdXNfYm9yZGVyIDwtIHRtX3NoYXBlKHVzX3N0YXRlc19zZikgKyB0bV9ib3JkZXJzKGNvbCA9ICJibGFjayIpICsgdG1fbGF5b3V0KGZyYW1lID0gRikgDQp1c19iZyArIHRvcm5hZG8gKyB1c19ib3JkZXINCmBgYA0KIyMjIENoZWNraW5nIENSUyBhbmQgZGF0YSBhdHRyaWJ1dGUNCmBgYHtyfQ0Kc3RfY3JzKHVzX3N0YXRlc19zZikNCmhlYWQodG9ybl9zZikNCnN1bW1hcnkodG9ybl9zZikNCnN0X2dlb21ldHJ5KHRvcm5fc2YpDQpgYGANCg0KIyMjIEF0dHJpYnV0ZSBzZWxlY3Rpb24NCmBgYHtyfQ0KaW5kZXggPC0gdXNfc3RhdGVzX3NmJFNUQVRFX05BTUUgPT0gIlRleGFzIiB8IHVzX3N0YXRlc19zZiRTVEFURV9OQU1FID09ICJOZXcgTWV4aWNvIiB8IA0KICAgICAgICAgdXNfc3RhdGVzX3NmJFNUQVRFX05BTUUgPT0gIk9rbGFob21hIiB8IHVzX3N0YXRlc19zZiRTVEFURV9OQU1FID09ICJBcmthbnNhcyIgIA0KDQpBb0lfc2YgPC0gdXNfc3RhdGVzX3NmW2luZGV4LF0NCkFvSV9iZyA8LSB0bV9zaGFwZShBb0lfc2YpICsgdG1fcG9seWdvbnMoImdyZXk5MCIpDQpBb0lfYmcgDQpgYGANCg0KIyMjIENsaXANCmBgYHtyfQ0KdG9ybl9jbGlwX3NmIDwtIHRvcm5fc2ZbQW9JX3NmLF0NCnRvcm5fY2xpcCA8LSB0bV9zaGFwZSh0b3JuX2NsaXBfc2YpICsgdG1fZG90cyhmaWxsID0gInJlZCIsIHNpemUgPSAwLjA0LCBmaWxsX2FscGhhID0gMC41KQ0KQW9JX2JnICsgdG9ybl9jbGlwDQoNCmhlYWQodG9ybl9jbGlwX3NmKQ0KYGBgDQoNCiMjIyBQb2ludC1Jbi1Qb2x5Z29uOiBVc2luZyBzdF9qb2luKCkNCmBgYHtyfQ0KQW9JX3Rvcm5fc2YgPC0gc3Rfam9pbih0b3JuX2NsaXBfc2YsIEFvSV9zZiwgam9pbiA9IHN0X3dpdGhpbikNCmhlYWQoQW9JX3Rvcm5fc2YpDQoNCkFvSV90b3JuX2RmIDwtIGFzLmRhdGEuZnJhbWUoQW9JX3Rvcm5fc2YpDQpoZWFkKEFvSV90b3JuX2RmKQ0Kc3VtbWFyeShBb0lfdG9ybl9kZikNCmBgYA0KDQojIyMgQ3JlYXRpbmcgbmV3IGRhdGEgZnJhbWUNCmBgYHtyfQ0KQW9JX3Rvcm5fZGYkU1RBVEVfTkFNRSA8LSBkcm9wbGV2ZWxzKEFvSV90b3JuX2RmJFNUQVRFX05BTUUpIA0KbmV3ZGY8LSBkYXRhLmZyYW1lKE5hbWUgPSBBb0lfdG9ybl9kZiRTVEFURV9OQU1FLERhbWFnZT0gQW9JX3Rvcm5fc2YkREFNQUdFKQ0KaGVhZChuZXdkZikNCm5ld2RmJERhbWFnZSA8LSBhcy5pbnRlZ2VyKG5ld2RmJERhbWFnZSkNCmBgYA0KDQojIyMgQ3Jvc3N0YWIgYW5hbHlzaXMNCmBgYHtyfQ0KY291bnQ8LSB0YWJsZShuZXdkZiREYW1hZ2UsIG5ld2RmJE5hbWUpDQpjb3VudA0KDQpjb3VudDI8LSB4dGFicyh+IG5ld2RmJERhbWFnZSArIG5ld2RmJE5hbSApDQpjb3VudDINCmBgYA0KDQojIyAxLjIgTWFwcGluZyBUb3JuYWRvIERlbnNpdHkgDQpgYGB7cn0NCiMgKENvdW50aW5nIFBvaW50cyBpbiBQb2x5Z29ucykNCg0KcG9pbnRzX3NmX2pvaW5lZCA8LSBzdF9qb2luKHRvcm5fc2YsIHVzX3N0YXRlc19zZiwgam9pbiA9IHN0X3dpdGhpbikNCmhlYWQocG9pbnRzX3NmX2pvaW5lZCkNCg0KdGFibGUxIDwtIHRhYmxlKHBvaW50c19zZl9qb2luZWQkU1RBVEVfTkFNRSkNCg0KY291bnRfc2YgPC0gYXMuZGF0YS5mcmFtZSh0YWJsZTEpDQpjb2xuYW1lcyhjb3VudF9zZikgPC0gYygiU1RBVEVfTkFNRSIsIkNvdW50cyIpDQpoZWFkKGNvdW50X3NmKQ0KDQpsaWJyYXJ5KGRwbHlyKQ0KdXNfc3RhdGVzX3NmPC0gbGVmdF9qb2luKHVzX3N0YXRlc19zZiwgY291bnRfc2YpDQpoZWFkKHVzX3N0YXRlc19zZikNCg0KYXJlYTE8LXN0X2FyZWEodXNfc3RhdGVzX3NmKSAjIHVuaXQ6IGZvb3QNCg0KbGlicmFyeSh1bml0cykNCmFyZWExPC1zZXRfdW5pdHMoYXJlYTEsIGttXjIpDQp1c19zdGF0ZXNfc2YkQVJFQTEgPC0gYXJlYTENCnVzX3N0YXRlc19zZiREZW5zaXR5IDwtIHVzX3N0YXRlc19zZiRDb3VudHMgLyB1c19zdGF0ZXNfc2YkQVJFQTENCmhlYWQodXNfc3RhdGVzX3NmKQ0KDQpgYGANCmBgYHtyfQ0KbGlicmFyeShwYWxzKSAjIGNvbGxlY3Rpb24gb2YgY29sb3IgcGFsZXR0ZXMNCnBsb3QodXNfc3RhdGVzX3NmWyJEZW5zaXR5Il0sIGJyZWFrcyA9ICJqZW5rcyIsIG5icmVha3MgPSA2LCBwYWw9YnJld2VyLnJlZHMoNikpDQoNCmBgYA0KDQojIEJ1ZmZlciB6b25lDQpgYGB7cn0NCnN0X2NycyhwbGFjZXNfc2YpDQpzdF9jcnMocGxhY2VzX3NmKTwtIHN0X2NycyhibG9ja3Nfc2YpDQpoZWFkKHBsYWNlc19zZikNCg0KYmxvY2tfYmcgPC0gdG1fc2hhcGUoYmxvY2tzX3NmKSArIHRtX3BvbHlnb25zKCJncmV5OTAiKQ0KcGxhY2VzX3B0cyA8LSB0bV9zaGFwZShwbGFjZXNfc2YpICsgdG1fZG90cyhmaWxsID0gInJlZCIsIHNpemUgPSAwLjUpDQoNCmJsb2NrX2JnICsgcGxhY2VzX3B0cw0KDQpwbGFjZXNfYnVmX3NmIDwtIHN0X2J1ZmZlcihwbGFjZXNfc2YsIGRpc3QgPSAzMDAwKQ0Kc3VtbWFyeShwbGFjZXNfYnVmX3NmKQ0KDQpwbGFjZXNfYnVmIDwtIHRtX3NoYXBlKHBsYWNlc19idWZfc2YpICsgdG1fcG9seWdvbnMoInllbGxvdyIpDQpibG9ja19iZyArIHBsYWNlc19idWYgKyBwbGFjZXNfcHRzDQoNCnBsYWNlc19idWZVX3NmIDwtIHN0X3VuaW9uKHBsYWNlc19idWZfc2YpDQpoZWFkKHBsYWNlc19idWZVX3NmKQ0KDQpwbGFjZXNfYnVmVSA8LSB0bV9zaGFwZShwbGFjZXNfYnVmVV9zZikgKyB0bV9wb2x5Z29ucygieWVsbG93IikNCmJsb2NrX2JnICsgcGxhY2VzX2J1ZlUgKyBwbGFjZXNfcHRzDQpgYGANCg0KIyBEaXN0YW5jZSBBbmFseXNpcw0KDQojIyBjaGVja2luZyBhdHRyaWJ1dGUgdGFibGUNCmBgYHtyfQ0KaGVhZChibG9ja3Nfc2YpDQpgYGANCg0KIyMgdXNpbmcgc3RfY2VudHJvaWQoKQ0KYGBge3J9DQpibG9ja3NfY2VudGVyX3NmIDwtIHN0X2NlbnRyb2lkKGJsb2Nrc19zZikNCmJsb2Nrc19jZW50ZXIgPC0gdG1fc2hhcGUoYmxvY2tzX2NlbnRlcl9zZikgKyB0bV9kb3RzKGZpbGwgPSAiYmx1ZSIsIHNpemUgPSAwLjIpDQpibG9ja19iZyArYmxvY2tzX2NlbnRlcg0KYmxvY2tzX2NlbnRlciArIHBsYWNlc19wdHMNCmBgYA0KDQojIyB1c2luZyBzdF9kaXN0YW5jZSgpDQpgYGB7cn0NCmRpc3RhbmNlLm1hdHJpeCA8LSBzdF9kaXN0YW5jZShibG9ja3NfY2VudGVyX3NmLCBwbGFjZXNfc2YpDQpsaWJyYXJ5KHVuaXRzKQ0KZGlzdGFuY2UubWF0cml4PC1zZXRfdW5pdHMoZGlzdGFuY2UubWF0cml4LCBrbSkNCmBgYA0KDQojIyB1c2luZyBhcHBseSgpDQpgYGB7cn0NCm5lYXJfZGlzdDwtIGFwcGx5KGRpc3RhbmNlLm1hdHJpeCwgMSwgbWVhbikNCm5lYXJfZGlzdFsxXQ0KDQp4aWQ8LSB3aGljaC5taW4obmVhcl9kaXN0KQ0KbmVhcl9kaXN0W3hpZF0NCmJsb2Nrc19zZlt4aWQsXQ0KYmxvY2tzX3NmJFBPUDE5OTBbeGlkXQ0KDQptaW5fZGlzdDwtIGFwcGx5KGRpc3RhbmNlLm1hdHJpeCwgMSwgbWluKQ0Kc2VsX2Jsb2NrczwtIG1pbl9kaXN0IDwgMSAgIyA8IDEga20NCnNlbF9zZiA8LSBibG9ja3Nfc2Zbc2VsX2Jsb2NrcyxdDQpzZWxfbWFwIDwtIHRtX3NoYXBlKHNlbF9zZikgKyB0bV9wb2x5Z29ucygieWVsbG93IikNCmJsb2NrX2JnICsgc2VsX21hcCArIHBsYWNlc19wdHMNCmBgYA0KDQojIyB1c2luZyBzdF9pc193aXRoaW5fZGlzdGFuY2UoKQ0KYGBge3J9DQpkMTwtIHNldF91bml0cygxLCBrbSkNCnNlbF9ibG9ja3MgPSBzdF9pc193aXRoaW5fZGlzdGFuY2UoYmxvY2tzX2NlbnRlcl9zZiwgcGxhY2VzX3NmLCBkaXN0ID0gZDEpDQpjbGFzcyhzZWxfYmxvY2tzKQ0Kc2VsX2Jsb2NrczwtIGxlbmd0aHMoc2VsX2Jsb2NrcykgPiAwDQpzZWxfc2YgPC0gYmxvY2tzX3NmW3NlbF9ibG9ja3MsXQ0Kc2VsX21hcCA8LSB0bV9zaGFwZShzZWxfc2YpICsgdG1fcG9seWdvbnMoInllbGxvdyIpDQpibG9ja19iZyArIHNlbF9tYXAgKyBwbGFjZXNfcHRzDQpgYGANCg0KIyMgY2hvb3NlIHNwZWNpZmljIG9uZSBwbGFjZQ0KYGBge3J9DQpzZWxfYmxvY2tzID0gc3RfaXNfd2l0aGluX2Rpc3RhbmNlKGJsb2Nrc19jZW50ZXJfc2YsIHBsYWNlc19zZlszLF0sIGRpc3QgPSBkMSkNCnNlbF9ibG9ja3MgPC0gbGVuZ3RocyhzZWxfYmxvY2tzKSA+IDANCnNlbF9zZiA8LSBibG9ja3Nfc2Zbc2VsX2Jsb2NrcyxdDQpzZWxfbWFwIDwtIHRtX3NoYXBlKHNlbF9zZikgKyB0bV9wb2x5Z29ucygieWVsbG93IikNCmJsb2NrX2JnICsgc2VsX21hcCArIHBsYWNlc19wdHMNCg0KVE9UQUxfUE9QIDwtIHN1bShzZWxfc2YkUE9QMTk5MCkNClRPVEFMX1dISVRFIDwtIHN1bShzZWxfc2YkUE9QMTk5MCAqIChzZWxfc2YkUF9XSElURS8xMDApICkNCmBgYA0KDQojIyBhY2Nlc3NpYmlsaXR5IGFuYWx5c2lzOiB1c2luZyBhcHBseSgpDQpgYGB7cn0NCmV0aDwtIGFzLmRhdGEuZnJhbWUoYmxvY2tzX2NlbnRlcl9zZlssMTQ6MThdKQ0KZXRoPC0gYXMubWF0cml4KGV0aFssMTo1XSkNCmV0aDwtIGFwcGx5KGV0aCwgMiwgZnVuY3Rpb24oeCkgKHgqYmxvY2tzX2NlbnRlcl9zZiRQT1AxOTkwKSkNCmNvbG5hbWVzKGV0aCk8LSBjKCJXaGl0ZSIsIkJsYWNrIiwiTmF0aXZlIiwiQXNpYW4iLCJPdGhlciIpDQpldGguYWNjZXNzIDwtIHh0YWJzKGV0aH5zZWxfYmxvY2tzKSANCg0KZXRoLmFjY2VzcyA8LSBhcy5kYXRhLmZyYW1lKGV0aC5hY2Nlc3MpDQpjb2xuYW1lcyhldGguYWNjZXNzKSA8LSBjKCJBY2Nlc3MiLCJFdGg1IiwiUG9wIikNCg0KbGlicmFyeShnZ3Bsb3QyKQ0KZ2dwbG90KGV0aC5hY2Nlc3MpICsgYWVzKGZpbGw9QWNjZXNzLCB5PVBvcCwgeD1FdGg1KSArIA0KICBnZW9tX2Jhcihwb3NpdGlvbj0ic3RhY2siLCBzdGF0PSJpZGVudGl0eSIpDQpgYGANCg0K