R新手求助：LiDAR CSV批量处理及数字场模型插值对比问题

Hey there! I totally get how overwhelming it can be when you’re new to R and tackling a big LiDAR project like this—let’s walk through each part step by step to get you sorted. You’ve got three main tasks: batch downloading your Litto3D CSV tiles, cleaning up column names and merging the data, then generating and comparing DFMs with different interpolation methods. Here’s how to do each one:

First, we’ll use the purrr package (part of the tidyverse) to loop through your tile URLs and download them all at once—way more efficient than manual downloads!

# Load tidyverse for easy looping and data handling
library(tidyverse)

# Replace this with your actual list of 12 Litto3D tile URLs
tile_urls <- c(
  "https://example.com/litto3d_tile_01.csv",
  "https://example.com/litto3d_tile_02.csv",
  # ... add the rest of your 10 tile URLs here
)

# Create a folder to store downloaded tiles (no error if it already exists)
output_folder <- "litto3d_lidar_tiles"
dir.create(output_folder, showWarnings = FALSE)

# Batch download each tile
walk(tile_urls, function(url) {
  # Extract filename from the URL
  tile_filename <- basename(url)
  # Download the file to your output folder
  download.file(url, file.path(output_folder, tile_filename), mode = "wb")
})

Litto3D CSVs often have non-standard column names (like X_LAMBER93, Z_ALTITUDE), so we’ll standardize them to X, Y, Z and merge all 12 tiles into one dataset.

# Get a list of all CSV files in your download folder
tile_files <- list.files(output_folder, pattern = "\\.csv$", full.names = TRUE)

# Read each tile, rename columns, and merge into a single dataframe
merged_lidar_data <- map_dfr(tile_files, function(file) {
  # Read the CSV—note: French datasets often use ; as separator and , as decimal
  tile_data <- read.csv(file, sep = ";", dec = ",")
  
  # Rename columns to standard X/Y/Z (replace with your actual column names!)
  tile_data <- tile_data %>%
    rename(
      X = X_LAMBER93,  # Swap with your actual X column name
      Y = Y_LAMBER93,  # Swap with your actual Y column name
      Z = Z_ALTITUDE   # Swap with your actual Z column name
    )
  
  # Keep only the columns we need for interpolation
  tile_data %>% select(X, Y, Z)
})

# Quick check to confirm the merge worked
head(merged_lidar_data)

We’ll use spatial packages like terra (modern replacement for raster) and gstat to run common interpolation methods, then visualize the differences.

# Load spatial packages
library(terra)
library(gstat)

# Convert merged dataframe to a spatial vector object
# Important: Litto3D uses EPSG:2154 (RGF93/Lambert-93) for French coordinates
spatial_lidar <- vect(merged_lidar_data, geom = c("X", "Y"), crs = "EPSG:2154")

# Create an empty raster template (adjust resolution to your needs—here it's 1m)
dfm_template <- rast(ext(spatial_lidar), resolution = 1, crs = crs(spatial_lidar))

## 3.1 Inverse Distance Weighting (IDW)
idw_dfm <- interpolate(dfm_template, spatial_lidar, method = "idw", formula = Z ~ 1)

## 3.2 Ordinary Kriging
# First, fit a variogram to model spatial autocorrelation
lidar_variogram <- variogram(Z ~ 1, spatial_lidar)
fitted_variogram <- fit.variogram(lidar_variogram, vgm("Sph"))  # Spherical model

# Run kriging interpolation
kriging_dfm <- interpolate(dfm_template, spatial_lidar, method = "krige", 
                           formula = Z ~ 1, model = fitted_variogram)

## 3.3 Nearest Neighbor Interpolation
nn_dfm <- interpolate(dfm_template, spatial_lidar, method = "near")

# Visualize all three DFMs side by side
par(mfrow = c(2, 2))
plot(idw_dfm, main = "IDW Interpolation")
plot(kriging_dfm, main = "Ordinary Kriging")
plot(nn_dfm, main = "Nearest Neighbor")
plot(spatial_lidar, add = TRUE, pch = 16, cex = 0.1)

Quick Tips for Your Project #

• Memory Management: LiDAR data is large! If you hit memory limits, try data.table instead of dplyr for faster, more memory-efficient reading/merging.
• CRS Double-Check: Always confirm your spatial data uses EPSG:2154 (Litto3D's standard) to avoid alignment errors.
• Interpolation Tweaks: For Kriging, test different variogram models (like vgm("Exp") for exponential). For IDW, adjust the power parameter with method = "idw", power = 2 (try 1 or 3 to see differences).

内容的提问来源于stack exchange，提问作者Ally U.