Spatial Analysis in Ecology

library(sf)

# point data
plot(st_point(c(0,0)))
plot(st_point(c(0.5,0.5)), add = T)

# line data
plot(st_linestring(matrix(c(0,0,1,1,2,0), ncol=2)) )

# polygon data
plot(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))))

library(sf)  # For handling spatial data
library(ggplot2)  # For plotting

# 五角星的polygon的顶点的数据
star_coords <- rbind(
  c(0, 1), c(0.2245, 0.309), c(0.9511, 0.309), c(0.3633, -0.118),
  c(0.5878, -0.809), c(0, -0.382), c(-0.5878, -0.809), c(-0.3633, -0.118),
  c(-0.9511, 0.309), c(-0.2245, 0.309), c(0, 1))

star_polygon <- st_polygon(list(star_coords))
star_sf <- st_sf(geometry = st_sfc(star_polygon, crs = 4326))

# 画图
plot(st_geometry(star_sf), col = "lightblue", border = "black", main = "Star Polygon")

> sh_polygon # 输入数据文件在R中的名字

Simple feature collection with 2 features and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 4 ymax: 4
CRS:           EPSG:4326
                        geometry
1 POLYGON ((0 0, 0 2, 2 2, 2 ...
2 POLYGON ((2 2, 2 4, 4 4, 4 ...

library(sf)  # For handling spatial data
library(ggplot2)  # For plotting
library(raster)  # For raster data manipulation

# Create a simple raster object
nrows <- 10  # Number of rows
ncols <- 10  # Number of columns
r <- raster(nrows = nrows, ncols = ncols, xmn = 0, xmx = 10, ymn = 0, ymx = 10)

# Assign random values to raster cells
set.seed(42)  # For reproducibility
values(r) <- runif(ncell(r), min = 0, max = 1)
raster_data <- as.data.frame(r, xy = TRUE)  # Convert raster to a data frame
raster_data[1:2,]  # Display the first two rows

> r # 栅格文件的名字

class      : RasterLayer 
dimensions : 10, 10, 100  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0.0002388966, 0.9888917  (min, max)

st_write(sf_object, "output.shp")
writeRaster(raster_object, "output.tif", format = "GTiff")

library(sf)
species_data <- data.frame(
  species = rep("Example_species", 10),
  # Random latitudes between -90 and 90
  latitude = runif(10, min = -90, max = 90),
  # Random longitudes between -180 and 180
  longitude = runif(10, min = -180, max = 180))

#write.csv(species_data, file = "species_points.csv",row.names=F) # save data
#species_data <- read.csv("species_points.csv") # read data

sf_data <- st_as_sf(species_data, coords = c("lon", "lat"), crs = 4326)
plot(sf_data)

library(sf)

# point shapefile (WGS84)
point <- st_as_sf(data.frame(lon = 120, lat = 30), 
                  coords = c("lon", "lat"), crs = 4326)
print(point)  # 查看经纬度坐标

# UTM projection (适合中国地区：带号50N)
point_utm <- st_transform(point, crs = 32650)
print(point_utm)  # 查看投影后的坐标

# crs = 32650 -- >
# crs = "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"

library(sf)  # For handling shapefiles and spatial data
library(ggplot2)  # For visualization
library(gridExtra)

# 1. 读取中国国家的 shapefile 数据
china <- st_read("spatial_data/CHN_adm0-3/CHN_adm0.shp", quiet = T)
st_crs(china) # 查看原始投影信息

# 2. 投影转换: 适合中国区域的投影 (Albers Equal Area)
albers_crs <- "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 
              +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
china_albers <- st_transform(china, crs = albers_crs)
st_crs(china_albers) # 查看转换后的投影信息

# 3. 可视化比较投影前后的数据
map_before <- ggplot() +
  geom_sf(data = china, fill = "lightblue", color = "black", alpha = 0.6) +
  ggtitle("原始投影") +
  theme_minimal()

map_after <- ggplot() +
  geom_sf(data = china_albers, fill = "lightgreen", color = "black", alpha = 0.6) +
  ggtitle("Albers Equal Area 投影") +
  theme_minimal()

grid.arrange(map_before, map_after)

library(sf)  # For handling shapefiles and spatial data
library(ggplot2)  # For visualization

# 3. 可视化比较投影前后的数据
map_before <- ggplot() +
  geom_sf(data = china, fill = "lightblue", color = "black", alpha = 0.6) +
  ggtitle("原始投影") +
  theme_minimal()

map_after <- ggplot() +
  geom_sf(data = china_albers, fill = "lightgreen", color = "black", alpha = 0.6) +
  ggtitle("Albers Equal Area 投影") +
  theme_minimal()

grid.arrange(map_before, map_after)

library(sf)
library(raster)  # For handling raster data

global_raster <- raster("bio01.tif") # 读取全球气候栅格文件
st_crs(global_raster) # 投影信息（通常为 WGS84）

# 定义目标投影 (Robinson)
robinson_crs <- "+proj=robin +lon_0=0 +datum=WGS84 +units=m +no_defs"

# 投影转换
global_raster_robinson <- projectRaster(global_raster, 
                      crs = robinson_crs, method = "bilinear")
crs(global_raster_robinson) # 查看转换后的投影信息

par(mfrow = c(1, 2))  # 将绘图窗口分为两部分
plot(global_raster, main = "原始投影 (WGS84)")
plot(global_raster_robinson, main = "Robinson 投影")

# 5. 保存转换后的栅格文件
writeRaster(global_raster_robinson, "global_climate_robinson.tif", format = "GTiff", overwrite = TRUE)

albers_crs <- "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 
              +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

robinson_crs <- "+proj=robin +lon_0=0 +datum=WGS84 +units=m +no_defs"

shp_subset <- shp_polygon %>% filter(field1 > 60)
shp_subset

st_write(shp_subset, "shp_subset", driver = "ESRI Shapefile", delete_layer = TRUE)
# 'delete_layer = TRUE' overwrite a shapefile with the same name

bio01 <- rast("spatial_data/bio01.tif") 
bio01
class(bio01)
names(bio01)
summary(bio01)
ext(bio01)
bio01_update <- bio01
bio01_update[ bio01_update < 0 ] = NA
plot(bio01_update)

bio01_update2 <- bio01 / 10
plot(bio01_update2)

library(terra)  # For raster and vector data handling
library(sf)     # For shapefile handling
bio01 <- rast("spatial_data/bio01.tif") 

# 根据随机点提取数据
ranSamps <- st_sample(st_as_sfc(st_bbox(bio01)), size=30, type = "random")
ranSamps <- st_as_sf(ranSamps) #Convert to simple feature collection
ranSamps
plot(bio01)
plot(ranSamps, add = T)

library(sf)
china_shp <- st_read("spatial_data/CHN_adm0-3/CHN_adm1.shp", quiet = TRUE)
head(as.data.frame(china_shp))
prov4_shp <- china_shp[china_shp$NAME_1 %in% c("Anhui","Jiangsu","Shanghai","Zhejiang"), ]
prov4_shp <- as(prov4_shp, "Spatial")
prov4_shp2 <- vect(prov4_shp)
prov4_shp2$value <- c("AH", "JS", "SH", "ZJ")
xy <- centroids(prov4_shp2)
plot(prov4_shp2, col=2:5)
points(xy, cex=4, pch=20, col='white')
text(prov4_shp2, 'value', cex=0.8)

library(sf)
china_shp <- st_read("spatial_data/CHN_adm0-3/CHN_adm1.shp", quiet = TRUE)
ginkgo_point <- read.csv("spatial_data/ginkgo_xy.csv")
ginkgo_point <- ginkgo_point %>% filter(complete.cases(.))
ginkgo_sf <- st_as_sf(ginkgo_point, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

# 裁剪
ginkgo_china <- st_crop(ginkgo_sf, china_shp)
plot(ginkgo_china)

library(sf)
china_shp <- st_read("spatial_data/CHN_adm0-3/CHN_adm1.shp", quiet = TRUE)
ginkgo_point <- read.csv("spatial_data/ginkgo_xy.csv")
ginkgo_point <- ginkgo_point %>% filter(complete.cases(.))
ginkgo_sf <- st_as_sf(ginkgo_point, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

#ginkgo_china <- st_crop(ginkgo_sf, china_shp)
#plot(ginkgo_china)
ginkgo_china <- st_intersection(ginkgo_sf, china_shp)

# need to make projection, convert to unit - m
albers_crs <- "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 
              +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
ginkgo_albers <- st_transform(ginkgo_china, crs = albers_crs)
# plot(st_geometry(ginkgo_albers))
china_albers <- st_transform(china_shp, crs = albers_crs)

# generate buffer around each point
buff_point <- st_buffer(ginkgo_albers, dist = 100000) # 10km-radius

ggplot() +
  geom_sf(data = china_albers, fill = "lightblue", color = "black") + 
  geom_sf(data = buff_point, fill = "darkgray", color = "black") + 
  geom_sf(data = ginkgo_albers, aes(color = "Points"), size = 1) + 
  theme_minimal() +
  theme(legend.position = "bottom")

# join the POINT geometries into a single MULTIPOINT geometry
ginkgo_union <- st_union(ginkgo_albers)
ginkgo_range <- st_convex_hull(ginkgo_union)
ginkgo_range_buff <- st_buffer(ginkgo_range, dist = 100000) # 10km-radius

ggplot() +
    geom_sf(data = ginkgo_range_buff, fill = NA, color = "black") + 
    geom_sf(data = ginkgo_range, fill = "lightblue", color = "lightblue") + 
    geom_sf(data = ginkgo_albers, aes(color = "Points"), size = 1) + 
    theme_minimal()

ginkgo_vect <- vect(ginkgo_albers)
# Define the extent and resolution of the raster
raster_extent <- ext(ginkgo_vect)  # Automatically match extent of points
resolution <- 50000  # Cell size (in map units, e.g., meters)

# Create an empty raster template
raster_template <- rast(ext = raster_extent, resolution = resolution, crs = crs(ginkgo_vect))

# Rasterize points: Count the number of points per cell
ginkgo_raster <- rasterize(ginkgo_vect, raster_template, fun = "length")

# Plot the raster
plot(ginkgo_raster, main = "Ginkgo Rasterized: Point Density")
plot(st_geometry(ginkgo_albers), add = TRUE, pch = 1, cex = 0.5)

r_df <- as.data.frame(ginkgo_raster, xy = TRUE)  # Extract cell values and coordinates

# Plot using ggplot2
ggplot() +
    # Add raster layer
    geom_tile(data = r_df, aes(x = x, y = y, fill = V1_length)) +
    scale_fill_viridis_c(name = "Value") +  # Color scale for raster values
    coord_equal() +  # Equal aspect ratio
    # Add spatial points
    geom_sf(data = ginkgo_albers, aes(geometry = geometry), 
            color = "red", size = 1, pch = 1) +
    # Add titles and themes
    labs(title = "Raster with Points Overlay",
         x = "Lon", y = "Lat") +
    theme_minimal()

library(terra)  # For raster and vector data handling
library(sf)     # For shapefile handling

bio01 <- rast("spatial_data/bio01.tif") 
china_shp <- st_read("spatial_data/CHN_adm0-3/CHN_adm0.shp", quiet = T)

# 需检查两类型的文件是否有相同投影系统
crs(bio01)
st_crs(china_shp)

# 裁剪/切割
cropped_raster <- crop(bio01, china_shp)  # 按范围裁剪
bio01_china <- mask(cropped_raster, vect(china_shp))  # 按多边形边界裁剪

par(mfrow = c(1, 2))  # 分为两个绘图区域
plot(bio01, main = "Original Raster")  # 原始栅格
plot(bio01_china, main = "Cropped Raster")  # 裁剪后的栅格

# 保存新的栅格文件
writeRaster(bio01_china, "bio01_cropped.tif", format = "GTiff", overwrite = TRUE)

# ref: https://rpubs.com/jonesor/gettinggbifdata
library(rgbif)
x <- name_backbone("species example")
gbif_download <- occ_download(
  pred("taxonKey", x$usageKey),
  pred("country", "DE"), 
  pred("hasGeospatialIssue", FALSE), 
  pred("hasCoordinate", TRUE), 
  pred("occurrenceStatus", "PRESENT"), 
  pred_gte("year", 1900), 
  user = user, pwd = pwd, email = email, # your account info
  format = "SIMPLE_CSV")

library(sf)
library(tidyverse)
library(rgbif)
#occ_data <- occ_search(taxonKey = name_backbone("Ginkgo biloba")$usageKey)
#head(occ_data$data)
#write.csv(occ_data$data, "spatial_data/ginkgo.csv", row.names = F, quote = F)
#ginkgo_point <- occ_data$data %>%
#    select(scientificName,decimalLatitude,decimalLongitude)
#write.csv(ginkgo_point, "spatial_data/ginkgo_xy.csv", row.names = F, quote = F)
ginkgo_point <- read.csv("spatial_data/ginkgo_xy.csv")
ginkgo_point <- ginkgo_point %>% filter(complete.cases(.))
ginkgo_sf <- st_as_sf(ginkgo_point, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)


library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

# 加载世界地图数据（中等分辨率）
world <- ne_countries(scale = "medium", returnclass = "sf")

ggplot() +
  geom_sf(data = world, fill = "lightblue", color = "black") +  # 世界地图
  geom_sf(data = ginkgo_sf, aes(color = "Points"), size = 1) +  # 点数据
  scale_color_manual(values = "red", labels = "Occurrence") +  # 自定义颜色
  theme_minimal() +
  ggtitle("World Map with Localities of Gingko") +
  theme(legend.position = "bottom")  # 调整图例位置

library(sf)
library(tidyverse)
library(rgbif)

# download data
occ_data <- occ_search(taxonKey = name_backbone("Ginkgo biloba")$usageKey)
head(occ_data$data)
#write.csv(occ_data$data, "spatial_data/ginkgo.csv", row.names = F, quote = F)
ginkgo_point <- occ_data$data %>%
    select(scientificName,decimalLatitude,decimalLongitude)
#write.csv(ginkgo_point, "spatial_data/ginkgo_xy.csv", row.names = F, quote = F)
ginkgo_point <- read.csv("spatial_data/ginkgo_xy.csv")
ginkgo_point <- ginkgo_point %>% filter(complete.cases(.))
ginkgo_sf <- st_as_sf(ginkgo_point, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

library(sf)
library(tidyverse)
library(rgbif)

# download data
# ...

library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

# 加载世界地图数据（中等分辨率）
world <- ne_countries(scale = "medium", returnclass = "sf")

ggplot() +
  geom_sf(data = world, fill = "lightblue", color = "black") +  # 世界地图
  geom_sf(data = ginkgo_sf, aes(color = "Points"), size = 1) +  # 点数据
  scale_color_manual(values = "red", labels = "Occurrence") +  # 自定义颜色
  theme_minimal() +
  ggtitle("World Map with Localities of Gingko") +
  theme(legend.position = "bottom")  # 调整图例位置

# extract bio01 values for points
bio01_values <- terra::extract(bio01_albers_cropped, ginkgo_albers)

ginkgo_sf_bio01 <- ginkgo_albers %>% mutate(values = bio01_values)
ginkgo_bio01_df <- st_coordinates(ginkgo_sf_bio01) %>%  # Extract x and y coordinates
    as.data.frame() %>%  # Convert to a data frame
    cbind(st_drop_geometry(ginkgo_sf_bio01))  # Bind other attributes (drop geometry column)

ginkgo_bio01_df2 <- st_coordinates(ginkgo_sf_bio01) %>%  # Extract x and y coordinates
    as.data.frame() %>%  # Convert to a data frame
    mutate(bio01 = bio01_values$bio01)

ginkgo_bio01_df2[1:3,]

library(spdep)  # For spatial autocorrelation
ginkgo_xy <- st_coordinates(ginkgo_albers)
ginkgo_bio01 <- bio01_values$bio01

# Create a neighbors list based on the coordinates
neighbors <- dnearneigh(ginkgo_xy, 0, 1000000)  # Specify a distance threshold (1000 km in projected CRS)
weights <- nb2listw(neighbors, style = "W")  # Create spatial weights

# Calculate Moran's I for bio01
moran_test <- moran.test(ginkgo_bio01, weights)

# Print Moran's I test result
print(moran_test)

# Visualization: Plot spatial autocorrelation
moran_plot <- moran.plot(ginkgo_bio01, listw = weights, main = "Moran's I Plot")