在R中基于SpatialPoint提取Landsat栅格值的问题与需求

Hey there! Let's fix up your workflow for extracting Landsat raster values to your points, plotting relationships, and exporting the data properly. I see where you're running into issues—let's break this down step by step.

First: Fix Coordinate Reference System (CRS) Mismatch #

A super common cause of extract() errors is mismatched CRS between your raster stack and SpatialPoints. Right now, your soil_sp doesn't have a CRS assigned, which will clash with your Landsat stack (which almost certainly has one set). Let's fix that first:

# Assign the same CRS from your Landsat stack to your SpatialPoints
proj4string(soil_sp) <- proj4string(lsat)

Correctly Extract and Merge Raster Values #

When you use extract() on a raster stack, it returns a matrix where each row is a point, and each column is a raster band. Trying to stuff this entire matrix into a single column (soil_sp$ref) is why you're getting weird results—you need to split it into separate columns and merge it with your point attributes.

Here's how to do it properly:

# Extract raster values (returns a matrix)
extracted_bands <- extract(lsat, soil_sp)

# Convert the matrix to a data frame and name columns to match your bands
extracted_df <- as.data.frame(extracted_bands)
colnames(extracted_df) <- names(lsat) # Uses the original band names from your stack

# Combine your original point attributes (from soil.clean2) with the extracted raster values
# This creates a SpatialPointsDataFrame with all your data
soil_spdf <- SpatialPointsDataFrame(
  coords = soil_sp,
  data = cbind(soil.clean2, extracted_df)
)

Plotting Relationships Between Extracted Values and Observations #

Now that all your data is in a single SpatialPointsDataFrame, plotting is straightforward. If you want to see relationships between all your raster bands and your observed variables, use pairs():

# Plot all pairwise relationships (adjust the variables to match your observed column names)
pairs(~ your_observation_column + b1 + b2 + b3 + b4 + b5 + b6_1 + b6_2 + b7, 
      data = soil_spdf,
      main = "Relationships Between Observations and Landsat Bands")

If you want a single scatter plot (e.g., Band 4 vs. your observation):

plot(your_observation_column ~ b4, 
     data = soil_spdf,
     xlab = "Landsat Band 4 Reflectance",
     ylab = "Your Observed Value",
     pch = 16, col = "steelblue")

Exporting the Final Data #

You can export either a regular DataFrame or keep it as a SpatialPointsDataFrame (for GIS use):

# Export as a CSV (regular DataFrame)
final_df <- as.data.frame(soil_spdf)
write.csv(final_df, "soil_with_landsat_values.csv", row.names = FALSE)

# Export as a Shapefile (SpatialPointsDataFrame)
library(rgdal)
writeOGR(
  obj = soil_spdf,
  dsn = "output_directory", # Folder to save the shapefile
  layer = "soil_landsat_data", # Name of the shapefile
  driver = "ESRI Shapefile",
  overwrite_layer = TRUE
)

Why Your Original Code Wasn't Working #

When you assigned soil_sp$ref <- extract(...), you were storing an entire matrix inside a single attribute column of your SpatialPoints. This makes plotting impossible (since R doesn't know how to handle a matrix in a single column for formula plots) and breaks the one-to-one match between points and their values. By splitting the matrix into individual columns and merging with your original data, you fix this entirely.

内容的提问来源于stack exchange，提问作者morteza