library(sf)
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
tracts <- st_read(
"C:/Users/noahs/Downloads/ES252/stockton.shp"
)
## Reading layer `stockton' from data source
## `C:\Users\noahs\Downloads\ES252\stockton.shp' using driver `ESRI Shapefile'
## Simple feature collection with 86 features and 38 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -121.4197 ymin: 37.86748 xmax: -121.1839 ymax: 38.05876
## Geodetic CRS: NAD83
names(tracts)
## [1] "tl_2020_06" "tl_2020_07" "tl_2020_08" "tl_2020_09" "tl_2020_10"
## [6] "tl_2020_11" "tl_2020_12" "tl_2020_13" "tl_2020_14" "tl_2020_15"
## [11] "tl_2020_16" "tl_2020_17" "Sheet1__Ma" "Sheet1__GE" "Sheet1__F2"
## [16] "Sheet1__Ob" "Shape_Leng" "Shape_Area" "LST_trend_" "LST_pvalue"
## [21] "NDVI_trend" "NDVI_pvalu" "NDBI_trend" "NDBI_pvalu" "Impervious"
## [26] "Impervio_1" "Mean_Albed" "GEOID" "LST_trend1" "LST_pval_1"
## [31] "NDVI_tre_1" "NDVI_pva_1" "NDBI_tre_1" "NDBI_pva_1" "Impervio_2"
## [36] "Impervio_3" "Mean_Alb_1" "OBJECTID_1" "geometry"
library(readxl)
stats <- read_excel(
"C:/Users/noahs/Downloads/Stockton_Tract_Trend_Data.xlsx"
)
names(stats)
## [1] "GEOID" "LST_trend_C_per_year" "LST_pvalue"
## [4] "NDVI_trend_per_year" "NDVI_pvalue" "NDBI_trend_per_year"
## [7] "NDBI_pvalue" "Impervious_mean" "Impervious_gt70_pct"
## [10] "Mean_Albedo_mean"
head(stats)
## # A tibble: 6 × 10
## GEOID LST_trend_C_per_year LST_pvalue NDVI_trend_per_year NDVI_pvalue
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 06077000101 0.138 0.000157 0.00199 2.03e- 8
## 2 06077000102 0.135 0.000111 0.00187 2.07e-10
## 3 06077000300 0.137 0.0000489 0.00103 2.13e- 4
## 4 06077000401 0.164 0.00000686 0.000477 1.36e- 1
## 5 06077000402 0.164 0.0000133 0.0000441 8.90e- 1
## 6 06077000500 0.160 0.0000145 0.00116 4.29e- 5
## # ℹ 5 more variables: NDBI_trend_per_year <dbl>, NDBI_pvalue <dbl>,
## # Impervious_mean <dbl>, Impervious_gt70_pct <dbl>, Mean_Albedo_mean <dbl>
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
tracts_joined <- left_join(
tracts,
stats,
by = "GEOID"
)
names(tracts_joined)
## [1] "tl_2020_06" "tl_2020_07" "tl_2020_08"
## [4] "tl_2020_09" "tl_2020_10" "tl_2020_11"
## [7] "tl_2020_12" "tl_2020_13" "tl_2020_14"
## [10] "tl_2020_15" "tl_2020_16" "tl_2020_17"
## [13] "Sheet1__Ma" "Sheet1__GE" "Sheet1__F2"
## [16] "Sheet1__Ob" "Shape_Leng" "Shape_Area"
## [19] "LST_trend_" "LST_pvalue.x" "NDVI_trend"
## [22] "NDVI_pvalu" "NDBI_trend" "NDBI_pvalu"
## [25] "Impervious" "Impervio_1" "Mean_Albed"
## [28] "GEOID" "LST_trend1" "LST_pval_1"
## [31] "NDVI_tre_1" "NDVI_pva_1" "NDBI_tre_1"
## [34] "NDBI_pva_1" "Impervio_2" "Impervio_3"
## [37] "Mean_Alb_1" "OBJECTID_1" "LST_trend_C_per_year"
## [40] "LST_pvalue.y" "NDVI_trend_per_year" "NDVI_pvalue"
## [43] "NDBI_trend_per_year" "NDBI_pvalue" "Impervious_mean"
## [46] "Impervious_gt70_pct" "Mean_Albedo_mean" "geometry"
summary(tracts_joined$LST_trend_C_per_year)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.009624 0.125617 0.143765 0.140637 0.165664 0.247559
tracts_joined <- st_transform(tracts_joined, 3310)
st_crs(tracts_joined)
## Coordinate Reference System:
## User input: EPSG:3310
## wkt:
## PROJCRS["NAD83 / California Albers",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["California Albers",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-120,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",34,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",40.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",-4000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["State-wide spatial data management."],
## AREA["United States (USA) - California."],
## BBOX[32.53,-124.45,42.01,-114.12]],
## ID["EPSG",3310]]
library(dplyr)
library(spdep)
## Warning: package 'spdep' was built under R version 4.5.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.5.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
tracts_moran <- tracts_joined %>%
filter(!is.na(LST_trend_C_per_year))
neighbors <- poly2nb(
tracts_moran,
queen = TRUE,
snap = 100
)
weights <- nb2listw(
neighbors,
style = "W",
zero.policy = TRUE
)
moran.test(
tracts_moran$LST_trend_C_per_year,
weights,
zero.policy = TRUE
)
##
## Moran I test under randomisation
##
## data: tracts_moran$LST_trend_C_per_year
## weights: weights
##
## Moran I statistic standard deviate = 4.1069, p-value = 2.005e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.256504471 -0.011764706 0.004266849
tracts_ndvi <- tracts_joined %>%
filter(!is.na(NDVI_trend_per_year))
neighbors_ndvi <- poly2nb(tracts_ndvi, queen = TRUE, snap = 100)
weights_ndvi <- nb2listw(
neighbors_ndvi,
style = "W",
zero.policy = TRUE
)
moran.test(
tracts_ndvi$NDVI_trend_per_year,
weights_ndvi,
zero.policy = TRUE
)
##
## Moran I test under randomisation
##
## data: tracts_ndvi$NDVI_trend_per_year
## weights: weights_ndvi
##
## Moran I statistic standard deviate = 2.3904, p-value = 0.008415
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.143510502 -0.011764706 0.004219564
# Remove NA values
tracts_ndbi <- tracts_joined %>%
filter(!is.na(NDBI_trend_per_year))
# Create spatial neighbors
neighbors_ndbi <- poly2nb(
tracts_ndbi,
queen = TRUE,
snap = 100
)
# Create spatial weights
weights_ndbi <- nb2listw(
neighbors_ndbi,
style = "W",
zero.policy = TRUE
)
# Run Global Moran's I
moran.test(
tracts_ndbi$NDBI_trend_per_year,
weights_ndbi,
zero.policy = TRUE
)
##
## Moran I test under randomisation
##
## data: tracts_ndbi$NDBI_trend_per_year
## weights: weights_ndbi
##
## Moran I statistic standard deviate = 5.0956, p-value = 1.739e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.319334124 -0.011764706 0.004222129