Spatial Data Science
Fardla
2025-10-22
Areal Data
the domain is a fixed countable collection of areal units at which variables are observed. For example, in spatial epidemiology, locations of individuals with a given desease are often aggregated in administrative areas. So it become like a sum or how much people in those areas by its administrative are infected of the desease.
# Data beras dan padi
padijateng <- read_excel("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/padijateng.xlsx")
# Batas wilayah kab/kota provinsi Jawa Tengah
indo <- st_read("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/gadm41_IDN_shp/gadm41_IDN_2.shp")## Reading layer `gadm41_IDN_2' from data source
## `C:\Users\WORKPLUS\Documents\Semester 5\Sains Data Spasial\gadm41_IDN_shp\gadm41_IDN_2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84jateng <- indo[indo$NAME_1 == "Jawa Tengah", ]
jateng <- jateng[-34,]
# Visualization
jateng_join <- jateng%>%left_join(padijateng, by = c("NAME_2" = "Kabupaten/Kota"))
ggplot(jateng_join) +
geom_sf(aes(fill = `Rekap Produksi Padi (ton)`), color = "white") +
scale_fill_viridis_c(option = "plasma", labels = label_number(big.mark = ".", decimal.mark = ",")) +
theme_minimal() +
labs(title = "Produksi Padi per Kabupaten/Kota di Jawa Tengah",
fill = "Produksi Padi (Ton)")
Geostatistical Data
D is a continuous fixed subset, so some location values can be predicted from the known spatial data and location. For example, air pollution measurement at a number of monitoring stations can be used to predict air pollution at other locations.
# Example for topsoil concentrations at several locations around the river Meuse
data(meuse)
meuse = st_as_sf(meuse, coords = c("x","y"), crs = 28992)
mapview(meuse, zcol = "lead", map.types = "CartoDB.Voyager")Point Pattern
The domain D is random. Its index set gives the locations of random events of the spatial point pattern, indicating occurance of the event or random giving some additional information. In this type of data the locations of corresponding event is even more important that the event itself.
# Data Gempa
datapoin <- read_excel("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/datagempa.xlsx")
names(datapoin) <- c("lintang","bujur","magnitudo")
# Buat datanya jadi geometri wilayah
datapoin_sf <- st_as_sf(datapoin,
coords = c("bujur", "lintang"),
crs = 4326)
# Peta indonesia
indo <- ne_countries(scale = "medium", returnclass = "sf") %>%
subset(admin == "Indonesia")
# Plot point pattern
ggplot() +
geom_sf(data = indo, fill = "antiquewhite", color = "#2C341B") +
geom_sf(data = datapoin_sf, color = "#FF006B", size = 2, alpha = 0.8) +
theme_minimal() +
labs(title = "Sebaran Titik Gempa di Indonesia",
subtitle = "Per Tanggal 28 Juli 2025",
x = "Bujur", y = "Lintang") +
theme(plot.title = element_text(size = 16, face = "bold"),
plot.subtitle = element_text(size = 12))
Vektor Data
Type of data that representing discreate geographical object using coordinate point
# Vektor Indonesia
ina <- st_read("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Provinsi Indonesia 38/Provinsi Indonesia 38/38_Provinsi_Indonesia.shp")## Reading layer `38_Provinsi_Indonesia' from data source
## `C:\Users\WORKPLUS\Documents\Semester 5\Sains Data Spasial\Provinsi Indonesia 38\Provinsi Indonesia 38\38_Provinsi_Indonesia.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 39 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range: zmin: -2.35e-05 zmax: -2.35e-05
## Geodetic CRS: WGS 84class(ina)## [1] "sf" "data.frame"print(ina)## Simple feature collection with 39 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range: zmin: -2.35e-05 zmax: -2.35e-05
## Geodetic CRS: WGS 84
## First 10 features:
## OBJECTID WADMPR Shape_Leng Shape_Area
## 1 1 <NA> 0.05864334 5.626854e-05
## 2 2 Aceh 29.25552992 4.630016e+00
## 3 3 Bali 6.10280874 4.589095e-01
## 4 4 Banten 10.87783440 7.645553e-01
## 5 5 Bengkulu 12.55618365 1.638437e+00
## 6 6 Daerah Istimewa Yogyakarta 3.46657637 2.599849e-01
## 7 7 DKI Jakarta 3.70407468 5.400618e-02
## 8 8 Gorontalo 11.28821860 9.770050e-01
## 9 9 Jambi 12.82777818 3.984796e+00
## 10 10 Jawa Barat 13.12718213 3.031137e+00
## geometry
## 1 MULTIPOLYGON (((98.17921 2....
## 2 MULTIPOLYGON (((97.38631 1....
## 3 MULTIPOLYGON (((115.1252 -8...
## 4 MULTIPOLYGON (((106.3499 -6...
## 5 MULTIPOLYGON (((102.385 -5....
## 6 MULTIPOLYGON (((110.8265 -8...
## 7 MULTIPOLYGON (((106.7934 -6...
## 8 MULTIPOLYGON (((121.3539 0....
## 9 MULTIPOLYGON (((104.4988 -1...
## 10 MULTIPOLYGON (((108.3231 -7...plot(ina)
pulau_jawa <- ina[c(4,6,7,10:12),]
plot(pulau_jawa)
Raster Data
Partioning the study areas to some polygon or rectangle area that have the same size. Where the partition is called pixel and on each pixel have a value of interest attribute
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
r## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev.tif
## name : elevation
## min value : 141
## max value : 547plot(r)
ncell(r)## [1] 8550Coordinate System
In the spatial analysis, location point is represented by coordinate and the coordinate itself is represented by the coordinate system. Generally, there’s two types of coordinate system, there are Geographic Coordinate System and Projection Coordinate System. GCS is build form the longitude and lattitude components, that is used for identificating location on a 3 dimensional earth, while the PCS is transform that to a 2 dimensional area in a meter;feet or others unit.
# Data vektor indonesia
ina <- st_make_valid(ina)
st_crs(ina)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## 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]]ina_albers <- st_transform(ina, crs = "+proj=aea +lat_1=-5 +lat_2=-15 +lat_0=0 +lon_0=120 +datum=WGS84 +units=m +no_defs")
st_crs(ina_albers)## Coordinate Reference System:
## User input: +proj=aea +lat_1=-5 +lat_2=-15 +lat_0=0 +lon_0=120 +datum=WGS84 +units=m +no_defs
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## 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",-5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",-15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]# Peta indonesia sistem koordinat geografis
tm_shape(ina) + tm_borders() + tm_graticules(lines = TRUE) + tm_scalebar(position = c("left","bottom")) + tm_compass(position = c("right","top"))
# Peta indonesia sistem koordinat proyeksi
tm_shape(ina_albers) + tm_borders(col = "blue") + tm_title("Koordinat Proyeksi (Albers Equal Area)") + tm_scalebar(position = c("left", "bottom")) + tm_compass(type = "8star", position = c("right", "top")) + tm_grid()
# Data raster
crs(r)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"# Transformasi sistem koordinat
r2 = terra::project(r, "EPSG:2169")
crs(r2)## [1] "PROJCRS[\"LUREF / Luxembourg TM\",\n BASEGEOGCRS[\"LUREF\",\n DATUM[\"Luxembourg Reference Frame\",\n ELLIPSOID[\"International 1924\",6378388,297,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4181]],\n CONVERSION[\"Luxembourg TM\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49.8333333333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",6.16666666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",1,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",80000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"northing (X)\",north,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"easting (Y)\",east,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Luxembourg.\"],\n BBOX[49.44,5.73,50.19,6.53]],\n ID[\"EPSG\",2169]]"Global Spatial Autocorrelation
Spatial autocorrelation is used to describe the extent to which a variable is correlated with itself through space. Positive spatial autocorrelation occurs when observations with similar values are closer together. Negative spatial autocorrelation occurs when observations with dissimilar values are closer together. Indices that usually used to assess spatial autocorrelation in areal data is Moran’s I. Moran’s I values usually range from –1 to 1. Moran’s I values significantly above E[I] indicate positive spatial autocorrelation or clustering. This occurs when neighboring regions tend to have similar values. Moran’s I values significantly below E[I] indicate negative spatial autocorrelation or dispersion. This happens when regions that are close to one another tend to have different values. Finally, Moran’s I values around E[I] indicate randomness, that is, absence of spatial pattern.
Local Moran’s I allows us to identify clusters of the following types: • High-High: areas of high values with neighbors of high values, • High-Low: areas of high values with neighbors of low values, • Low-High: areas of low values with neighbors of high values, • Low-Low: areas of low values with neighbors of low values.
data<-read_excel('C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Dataset-20250915/Data Jawa Timur.xlsx')
str(data)## tibble [38 × 6] (S3: tbl_df/tbl/data.frame)
## $ Kab/Kota : chr [1:38] "Pacitan" "Ponorogo" "Trenggalek" "Tulungagung" ...
## $ lintang : num [1:38] -8.1 -7.97 -8.14 -8.08 -8.14 -7.8 -8.09 -8.12 -8.26 -8.33 ...
## $ bujur : num [1:38] 111 112 112 112 112 ...
## $ Pert.Kendaraan: num [1:38] 189426 508956 322612 743627 647213 ...
## $ IPM2020 : num [1:38] 68.4 70.8 69.7 73 70.6 ...
## $ DBD : num [1:38] 84 42 144 129 108 32 152 65 59 31 ...# file shp
jatim <- st_read("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/File shp -20250915/poly3500.shp")## Reading layer `poly3500' from data source
## `C:\Users\WORKPLUS\Documents\Semester 5\Sains Data Spasial\File shp -20250915\poly3500.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 110.8987 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## CRS: NA# input data DBD ke dalam peta
DBD = data[1:38,6]
DBD.matrix = as.matrix(DBD)
data.gab = data.frame(jatim,DBD.matrix)
sd = merge(jatim,data.gab)
# ubah ke objek spasial
jatim_sp <- as(sd, "Spatial")
spplot(jatim_sp, colnames(DBD.matrix),scales=list(draw=F),colorkey=list(space="bottom"))
# Uji autokorelasi spasial global data DBD
nb <- poly2nb(jatim_sp, queen = TRUE)## Warning in poly2nb(jatim_sp, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style = "W")
moran.test(jatim_sp@data$DBD, lw)##
## Moran I test under randomisation
##
## data: jatim_sp@data$DBD
## weights: lw
##
## Moran I statistic standard deviate = 2.1038, p-value = 0.0177
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.24127269 -0.02702703 0.01626419# input data IPM ke dalam peta
IPM<-data[1:38,5]
IPM.matrix = as.matrix(IPM)
data.gab<-data.frame(jatim,IPM.matrix)
sd<-merge(jatim,data.gab)
# ubah ke objek spasial
jatim_sp <- as(sd, "Spatial")
spplot(jatim_sp, colnames(IPM.matrix),scales=list(draw=F),colorkey=list(space="bottom"))
# Uji autokorelasi spasial global data IPM
nb <- poly2nb(jatim_sp, queen = TRUE)## Warning in poly2nb(jatim_sp, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style = "W")
moran.test(jatim_sp@data$IPM, lw)##
## Moran I test under randomisation
##
## data: jatim_sp@data$IPM
## weights: lw
##
## Moran I statistic standard deviate = 3.4178, p-value = 0.0003156
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.42072876 -0.02702703 0.01716247Local Spatial Autocorrelation
There is often interest in providing a local measure of similarity between each area’s value and those of nearby areas. Local Indicators of Spatial Association (LISA) are designed to provide an indication of the extent of significant spatial clustering of similar values around each observation. Typically, the values of the LISA are mapped to indicate the location of areas with comparatively high or low local association with neighboring areas. A high value for Ii suggests that the area is surrounded by areas with similar values. Such an area is part of a cluster of high observations, low observations, or moderate observations. A low value for Ii indicates that the area is surrounded by areas with dissimilar values. Such an area is an outlier indicating that the observation of area i is different from most or all of the observations of its neighbors.
# file shp
jatim <- st_read("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/File shp -20250915/poly3500.shp")## Reading layer `poly3500' from data source
## `C:\Users\WORKPLUS\Documents\Semester 5\Sains Data Spasial\File shp -20250915\poly3500.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 110.8987 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## CRS: NA# dataset
data<-read_excel('C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Dataset-20250915/Data Jawa Timur.xlsx')
IPM <- data[1:38,5]
IPM.matrix = as.matrix(IPM)
data.gab <- data.frame(jatim,IPM.matrix)
jatim_sf <- merge(jatim,data.gab)
tm_shape(jatim_sf) + tm_fill("IPM", fill.scale = tm_scale_continuous(values = c("yellow","white","green"),
n =9),
fill.legend = tm_legend(title = "IPM")) + tm_borders(col = "black",lwd = 0.5) + tm_title("IPM Provinsi Jawa Timur 2020") + tm_layout(legend.outside = TRUE)## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
nb <- poly2nb(jatim_sf)## Warning in poly2nb(jatim_sf): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style = "W")
mp <- moran.plot(as.vector(scale(jatim_sf$IPM)), lw)
# Moran's lokal
loc_moran <- localmoran(jatim_sf$IPM, lw, alternative = "greater")
head(loc_moran)## Ii E.Ii Var.Ii Z.Ii Pr(z > E(Ii))
## 1 0.009821889 -8.630909e-04 0.0047029795 0.1558071 0.4380925
## 2 0.207549235 -1.670357e-02 0.3033981183 0.4071284 0.3419568
## 3 0.091456892 -3.297814e-03 0.0393215233 0.4778434 0.3163808
## 4 -0.021631491 -8.697760e-04 0.0058707092 -0.2709680 0.6067922
## 5 -0.044557417 -9.182763e-04 0.0079893131 -0.4882266 0.6873053
## 6 -0.006321085 -3.920977e-05 0.0002138301 -0.4295907 0.6662533jatim_sf$lmI <- loc_moran[, "Ii"]
jatim_sf$lmZ <- loc_moran[,"Z.Ii"]
jatim_sf$lmp <- loc_moran[,"Pr(z > E(Ii))"]
tmap_mode("plot")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`p1 <- tm_shape(jatim_sf) + tm_polygons(fill = "IPM", fill.scale = tm_scale_intervals(style = "quantile",
values = c("yellow","white","green")),
fill.legend = tm_legend(title = "IPM")) + tm_layout(legend.outside = TRUE)
p1
p2 <- tm_shape(jatim_sf) +
tm_polygons(fill = "lmI", fill.scale = tm_scale_intervals(style = "quantile", values = "brewer.blues",
midpoint = NA),
fill.legend = tm_legend(title = "Local Moran's I")) +
tm_layout(legend.outside = TRUE)
p2
p3 <- tm_shape(jatim_sf) +
tm_polygons(fill = "lmZ", fill.scale = tm_scale_intervals(breaks = c(-Inf,1.65,Inf),labels = c("Tidak Signifikan","Signifikan"),
values = c("grey80","orange"),
midpoint = NA),fill.legend = tm_legend(title = "Z-score")) +tm_layout(legend.outside = TRUE)
p3
p4 <- tm_shape(jatim_sf) +
tm_polygons(fill = "lmp", fill.scale = tm_scale_intervals(breaks = c(-Inf,0.05,Inf),labels = c("Signifikan","Tidak Signifikan"),
values = c("purple","grey80"),
midpoint = NA),
fill.legend = tm_legend(title = "p-value")) +
tm_layout(legend.outside = TRUE)
p4
tmap_arrange(p1,p2,p3,p4)
jatim_sf$quadrant <- NA
jatim_sf[(mp$x >= 0 & mp$wx >= 0) & (jatim_sf$lmp <= 0.05), "quadrant"] <- 1
jatim_sf[(mp$x <= 0 & mp$wx <= 0) & (jatim_sf$lmp <= 0.05), "quadrant"] <- 2
jatim_sf[(mp$x >= 0 & mp$wx <= 0) & (jatim_sf$lmp <= 0.05), "quadrant"] <- 3
jatim_sf[(mp$x <= 0 & mp$wx >= 0) & (jatim_sf$lmp <= 0.05), "quadrant"] <- 4
jatim_sf[(jatim_sf$lmp > 0.05),"quadrant"] <- 5
p_quadrant <- tm_shape(jatim_sf) + tm_polygons(fill = "quadrant",
fill.scale = tm_scale_intervals(breaks = c(1,2,3,4,5,6),
labels = c("High-High","Low-Low","High-Low","Low-High", "Non-significant"),
values = c("green","red","yellow","orange","white")),
fill.legend = tm_legend(title = "")) + tm_borders(fill_alpha = 0.5) + tm_title("Clusters") + tm_layout(frame = FALSE, legend.outside = TRUE)
p_quadrant## Multiple palettes called "green" found: "kovesi.green", "tableau.green". The first one, "kovesi.green", is returned.
Spatial Data Operation
# File batas administratif desa di Indonesia
indo_desa <- st_read("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/BATAS WILAYAH KELURAHAN-DESA 10K/Batas_Wilayah_KelurahanDesa_10K_AR.shp")## Reading layer `Batas_Wilayah_KelurahanDesa_10K_AR' from data source
## `C:\Users\WORKPLUS\Documents\Semester 5\Sains Data Spasial\BATAS WILAYAH KELURAHAN-DESA 10K\Batas_Wilayah_KelurahanDesa_10K_AR.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 83518 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range: zmin: 0 zmax: 1876.818
## Geodetic CRS: WGS 84names(indo_desa)## [1] "OBJECTID" "NAMOBJ" "FCODE" "REMARK" "METADATA" "SRS_ID"
## [7] "KDBBPS" "KDCBPS" "KDCPUM" "KDEBPS" "KDEPUM" "KDPBPS"
## [13] "KDPKAB" "KDPPUM" "LUASWH" "TIPADM" "WADMKC" "WADMKD"
## [19] "WADMKK" "WADMPR" "WIADKC" "WIADKK" "WIADPR" "WIADKD"
## [25] "UUPP" "LUAS" "geometry"unique(indo_desa$WADMPR)## [1] "Jawa Barat" "Bali"
## [3] "Nusa Tenggara Barat" "Jawa Tengah"
## [5] "Daerah Istimewa Yogyakarta" "Jawa Timur"
## [7] "Sulawesi Selatan" "Sumatera Barat"
## [9] "Jambi" "Lampung"
## [11] "Sumatera Selatan" "Bengkulu"
## [13] "Banten" "Kepulauan Bangka Belitung"
## [15] "DKI Jakarta" "Kalimantan Selatan"
## [17] "Kalimantan Barat" "Kalimantan Tengah"
## [19] "Nusa Tenggara Timur" "Maluku"
## [21] "Papua Selatan" "Sulawesi Barat"
## [23] "Papua Tengah" "Papua Pegunungan"
## [25] "Papua Barat" "Papua"
## [27] "Sulawesi Tengah" "Sulawesi Tenggara"
## [29] "Kepulauan Riau" "Riau"
## [31] "Sumatera Utara" "Kalimantan Timur"
## [33] "Kalimantan Utara" "Sulawesi Utara"
## [35] "Aceh" NA
## [37] "Maluku Utara" "Papua Barat Daya"
## [39] "Gorontalo"# Subset Jawa Tengah
jateng = indo_desa %>% filter(WADMPR == "Jawa Tengah")
plot(st_geometry(jateng), main = "Wilayah Jawa Tengah (batas desa)")
# Ubah geometri dari 3D ke 2D
jateng <- st_zm(jateng, drop = TRUE, what = "ZM")
# Pastikan CRS tetap ada
st_crs(jateng) <- 4326
# Menyimpan shapefile
dir.create("C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Peta Jateng1")
st_write(jateng, "C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Peta Jateng1/Batas Wilayah Kelurahan-Desa Jateng.shp")## Writing layer `Batas Wilayah Kelurahan-Desa Jateng' to data source
## `C:/Users/WORKPLUS/Documents/Semester 5/Sains Data Spasial/Peta Jateng1/Batas Wilayah Kelurahan-Desa Jateng.shp' using driver `ESRI Shapefile'
## Writing 8562 features with 26 fields and geometry type Multi Polygon.# Cek class data Jateng
class(jateng)## [1] "sf" "data.frame"st_crs(jateng)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]jateng_proj <- st_transform(jateng, 32749)
# Menggabungkan desa menjadi kabupaten dan kota
jateng_kab.kota <- jateng_proj %>% group_by(WADMKK) %>% summarise(geometry = st_union(geometry))
plot(st_geometry(jateng_kab.kota), col = "lightblue", main = "Union (Dissolve) Jawa Tengah")
# Membuat centroid dari tiap kabupaten/kota
jateng_centroid <- st_centroid(jateng_kab.kota)## Warning: st_centroid assumes attributes are constant over geometries# Menampilkan centroid di atas peta kabupaten/kota
plot(st_geometry(jateng_kab.kota),
col = "lightblue",
main = "Centroid Setiap Kabupaten/Kota di Jawa Tengah")
# Tambahkan titik centroid (warna biru)
plot(st_geometry(jateng_centroid),
col = "red",
pch = 19,
add = TRUE)
# Data titik bandara (maps)
bandaradf <- data.frame(
Nama = c("Ahmad Yani", "Adi Soemarno","Tunggul Wulung"),
lon = c(110.374,110.753,109.032),
lat = c(-6.971,-7.514,-7.644)
)
# Konversi ke data spasial
bandarasf <- st_as_sf(bandaradf, coords = c("lon","lat"), crs = 4326)
# Memastikan proyeksi sudah sama
jateng_kab.kota = st_transform(jateng_kab.kota, 32749)
bandarasf = st_transform(bandarasf, 32749)
# Membuat buffer 30 km di sekitar bandara
bandarabuffer = st_buffer(bandarasf, dist = 30000)
ggplot() + geom_sf(data = jateng_kab.kota, fill = "white", color = "grey50") + geom_sf(data = bandarabuffer, fill = "lightblue", alpha =
0.5) + geom_sf(data = bandarasf, color = "red", size = 3) + geom_sf_text(data = bandarasf, aes(label = Nama), nudge_y = 20000, size = 3.5) + ggtitle("Zona Buffer 30 km di Sekitar Bandara Jawa Tengah") + theme_minimal() 
# Overlay kabupaten yang masuk dalam zona buffer bandara (intersection)
kab_dalam_buffer <- st_intersection(jateng_kab.kota, bandarabuffer)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesggplot() + geom_sf(data = jateng_kab.kota, fill = "white", color = "grey70") + geom_sf(data = kab_dalam_buffer, aes(fill = Nama), alpha = 0.5) + geom_sf(data = bandarabuffer, color = "blue", fill = NA, linetype = "dashed") + geom_sf(data = bandarasf, color = "red", size = 3) + geom_sf_text(data = bandarasf, aes(label = Nama), nudge_y = 20000) + ggtitle("Kabupaten yang Masuk dalam Zona Buffer 30 km Bandara di Jawa Tengah") + scale_fill_viridis_d(option = "C", name = "Kabupaten") + theme_minimal()
# Menghitung luas wilayah kab/kota yang masuk dalam zona buffer
kab_dalam_buffer$luas_km2 <- as.numeric(st_area(kab_dalam_buffer)) / 1e6
hasil_ringkasan <- kab_dalam_buffer %>% st_drop_geometry() %>% group_by(WADMKK) %>% summarise(total_luas_km2 = sum(luas_km2, na.rm = TRUE)) %>% arrange(desc(total_luas_km2))
hasil_ringkasan## # A tibble: 16 × 2
## WADMKK total_luas_km2
## <chr> <dbl>
## 1 Cilacap 1222.
## 2 Boyolali 883.
## 3 Kendal 669.
## 4 Banyumas 663.
## 5 Semarang 533.
## 6 Sragen 513.
## 7 Klaten 467.
## 8 Demak 438.
## 9 Sukoharjo 430.
## 10 Kota Semarang 370.
## 11 Karanganyar 286.
## 12 Grobogan 64.3
## 13 Kota Surakarta 46.7
## 14 Temanggung 25.2
## 15 Wonogiri 1.70
## 16 Kota Salatiga 0.543# Operasi CROP / CLIP
jateng_clip <- st_intersection(jateng_kab.kota, bandarabuffer)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Menambahkan kolom luas hasil potongan (opsional)
jateng_clip$luas_km2 <- as.numeric(st_area(jateng_clip)) / 1e6
# --- Visualisasi hasilnya ---
ggplot() + geom_sf(data = jateng_kab.kota, fill = "white", color = "grey70") + geom_sf(data = bandarabuffer, fill = "lightblue", alpha = 0.4) + geom_sf(data = jateng_clip, fill = "tomato", alpha = 0.6, color = "black") + ggtitle("Hasil Crop (Clip) Wilayah Jawa Tengah dengan Zona Buffer Bandara") + theme_minimal()