感谢你提供的列名信息，我们可以看到：
- `"...1"`：可能是 Excel 中第一列无标题的自动命名（比如行号或空列），可以忽略。
- `"STATION"`：站点编号
- `"NAME"`：站点名称
- `"LATITUDE"`：纬度
- `"LONGITUDE"`：经度
- `"ELEVATION"`：海拔
- `"statement"`：可能是一个状态或注释字段
- 后续从 `"2020-01-01"` 到 `"2020-01-31"` 是每日风速数据（共31天）
---
### ✅ 目标更新：
你要做的是：**根据给定的 SHP 范围裁剪站点，保留落在该地理范围内的所有站点及其完整的逐日风速数据。**
下面是适配你实际列名的完整 R 语言代码，并处理好坐标、投影和数据结构问题。
```r
# 加载所需包
library(sf)
library(readxl)
library(dplyr)
# ------------------- 参数设置 -------------------
excel_file <- "your_wind_station_data.xlsx" # 替换为你的实际文件路径
sheet_name <- "Sheet1" # 替换为你的工作表名
shp_file <- "your_boundary.shp" # 替换为你的SHP文件路径
# ------------------- 步骤1: 读取Excel数据 -------------------
df <- read_excel(excel_file, sheet = sheet_name)
# 查看列名确认
names(df) <- make.names(names(df)) # 确保列名是合法的（防止空格等问题）
cat("原始列名:\n"); print(names(df))
# ------------------- 步骤2: 提取空间信息并创建sf对象 -------------------
# 使用 LATITUDE 和 LONGITUDE 创建空间点（注意：顺序是 LON, LAT）
stations_sf <- st_as_sf(df,
coords = c("LONGITUDE", "LATITUDE"), # 经度在前，纬度在后
crs = 4326, # WGS84 地理坐标系
dim = "XY")
# ------------------- 步骤3: 读取SHP边界并确保其CRS -------------------
boundary <- st_read(shp_file)
# 如果boundary不是投影坐标系（如EPSG:32649），则需要检查并转换
# 假设你知道目标投影是 UTM Zone 49N (EPSG:32649)，我们统一到这个坐标系进行空间操作
if (is.na(st_crs(boundary))) {
stop("SHP文件没有坐标系信息，请先定义正确的CRS！")
}
# 将站点数据重投影到与SHP相同的坐标系下进行空间判断
stations_projected <- st_transform(stations_sf, crs = st_crs(boundary))
# ------------------- 步骤4: 空间裁剪 —— 找出在SHP范围内的站点 -------------------
# 使用 st_intersects 或 st_within 进行空间子集提取
# 这里使用 [s,t] 语法：返回落在任意多边形内的站点
stations_clipped <- stations_projected[boundary, , op = st_intersects]
# 若你想更严格地要求“完全包含”，可用 st_within，但通常 st_intersects 更通用
# ------------------- 步骤5: 转回原始经纬度并提取属性表格 -------------------
# 将结果转回WGS84以便保留原始经纬度格式输出
stations_wgs84 <- st_transform(stations_clipped, crs = 4326)
# 去掉geometry列，恢复为普通data.frame，同时保留原始所有列（包括每日数据）
result_df <- st_drop_geometry(stations_wgs84)
# 可选：重新排序列，把时间序列放在后面
date_cols <- grep("^\\d{4}-\\d{2}-\\d{2}", names(result_df), value = TRUE)
non_date_cols <- setdiff(names(result_df), date_cols)
final_df <- select(result_df, c(non_date_cols, date_cols)) # 按逻辑排序
# ------------------- 步骤6: 导出结果 -------------------
write.csv(final_df, "clipped_stations_202001.csv", row.names = FALSE, na = "")
cat("共保留了", nrow(final_df), "个站点在SHP范围内。\n")
```
---
### ✅ 关键说明：
| 功能 | 说明 |
|------|------|
| `coords = c("LONGITUDE", "LATITUDE")` | 必须是 **经度在前，纬度在后**，否则位置错误！ |
| `crs = 4326` | 表示输入的经纬度使用 WGS84 坐标系 |
| `st_transform(...)` | 将点从地理坐标（度）转为投影坐标（米），确保与 SHP 在同一空间参考下比较 |
| `st_intersects` | 判断点是否与多边形相交（即落在内部），适用于大多数情况 |
| `st_drop_geometry()` | 移除空间结构，得到纯数据框用于导出 |
> 💡 输出的 CSV 文件将包含原始所有列，包括 `STATION`, `NAME`, `LATITUDE`, `LONGITUDE`, `ELEVATION`, `statement` 和所有日期列（如 `2020-01-01` 等），仅保留位于 SHP 范围内的站点。
---
### ✅ 示例输出片段（final_df 头几行）：
```
...1 STATION NAME LATITUDE LONGITUDE ELEVATION statement 2020-01-01 2020-01-02 ...
1 1 101 Beijing 39.90 116.4 50.0 good 3.2 4.1
2 2 102 Tianjin 39.08 117.2 10.0 good 5.0 3.8
...
```
---
###
扫一扫
4449
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
