在Bing地图叠加eddy-covariance system 2D足迹结果的技术问询

Got it, let's walk through how to successfully overlay your 2D eddy-covariance footprint (plotted with image.plot) onto the Bing map you're creating. The key here is handling coordinate system conversions properly, since your footprint uses local meter-based coordinates while the map uses geographic (lat/lon) or Mercator projection.

Step 1: Load Required Packages #

First, make sure you have all the necessary libraries installed and loaded:

library(OpenStreetMap)
library(rgdal)
library(fields)  # For image.plot
library(raster)  # For easier raster handling (optional but helpful)

Step 2: Load Footprint Data & Define Grid Parameters #

Load your footprint data and set up the local grid coordinates based on your domain:

# Load footprint data
footprint_data <- readRDS("FFP.rds")

# Define domain and grid properties
domain <- c(-100, 1000, -100, 1000)
cell_size <- 2
x_local <- seq(domain[1], domain[2], by = cell_size)
y_local <- seq(domain[3], domain[4], by = cell_size)

Step 3: Convert Local Coordinates to Geographic (Lat/Lon) #

Your footprint is centered at (0,0) which corresponds to your site's lat/lon. We'll convert the local meter grid to lat/lon using UTM projection (since UTM uses meters, matching your local coordinates):

# Site geographic coordinates
site_lon <- -97.191391
site_lat <- 36.055935

# Define coordinate reference systems (CRS)
wgs84 <- CRS("+proj=longlat +datum=WGS84")  # Standard lat/lon
utm_zone <- UTMzone(site_lon, site_lat)  # Auto-get correct UTM zone (should be 15N here)
utm_crs <- CRS(paste0("+proj=utm +zone=", utm_zone, " +datum=WGS84 +units=m"))

# Convert site lat/lon to UTM coordinates
site_utm <- spTransform(
  SpatialPoints(data.frame(x = site_lon, y = site_lat), proj4string = wgs84),
  utm_crs
)

# Calculate UTM coordinates for every grid point
grid_utm_x <- site_utm@coords[1, 1] + x_local
grid_utm_y <- site_utm@coords[1, 2] + y_local

# Convert UTM grid back to lat/lon
grid_lonlat <- spTransform(
  SpatialPoints(expand.grid(grid_utm_x, grid_utm_y), proj4string = utm_crs),
  wgs84
)

# Reshape lat/lon into matrices matching your footprint data
grid_lon <- matrix(grid_lonlat@coords[, 1], nrow = length(x_local), ncol = length(y_local))
grid_lat <- matrix(grid_lonlat@coords[, 2], nrow = length(x_local), ncol = length(y_local))

Step 4: Fetch & Plot the Bing Map #

Fetch the map and convert it to WGS84 to match your lat/lon grid:

# Fetch Bing map (adjust bounds as needed)
map <- openmap(
  upperLeft = c(36.05778, -97.19250),
  lowerRight = c(36.05444, -97.18861),
  type = "bing"
)

# Convert map to WGS84 projection
map_wgs84 <- openproj(map, projection = wgs84)

# Plot the base map
plot(map_wgs84, main = "Eddy-Covariance Footprint Overlay")

Step 5: Overlay the Footprint Layer #

Now overlay your footprint with image and add a legend with image.plot:

# Plot footprint with transparency (adjust alpha for visibility)
image(
  x = grid_lon,
  y = grid_lat,
  z = footprint_data,
  add = TRUE,
  col = heat.colors(100),
  alpha = 0.6,  # Transparency to keep map visible underneath
  border = NA
)

# Add color legend (adjust position to avoid blocking map content)
image.plot(
  z = footprint_data,
  legend.only = TRUE,
  col = heat.colors(100),
  legend.width = 1,
  legend.shrink = 0.8,
  vertical = TRUE,
  x = par("usr")[2] + 0.02,  # Place legend to the right of the map
  y = mean(par("usr")[3:4])
)

Alternative: Using Raster Package for Easier Alignment #

If you prefer a more streamlined workflow, use the raster package to handle projection alignment automatically:

# Create a raster from your footprint data
footprint_raster <- raster(
  footprint_data,
  xmn = min(grid_lon),
  xmx = max(grid_lon),
  ymn = min(grid_lat),
  ymx = max(grid_lat),
  crs = wgs84
)

# Project raster to match the map's Mercator projection
footprint_merc <- projectRaster(footprint_raster, crs = projection(map_wgs84))

# Overlay on map
plot(footprint_merc, add = TRUE, col = heat.colors(100), alpha = 0.6, legend = FALSE)

# Add legend (same as above)
image.plot(
  z = footprint_data,
  legend.only = TRUE,
  col = heat.colors(100),
  legend.width = 1,
  legend.shrink = 0.8,
  vertical = TRUE,
  x = par("usr")[2] + 0.02,
  y = mean(par("usr")[3:4])
)

Key Notes to Troubleshoot #

• Matrix Dimension Match: Ensure footprint_data has dimensions matching length(x_local) × length(y_local). If not, transpose the matrix with t(footprint_data).
• UTM Zone Accuracy: Using UTMzone() ensures you get the correct zone for your site—this avoids coordinate shifts.
• Transparency: Adjust the alpha parameter to balance footprint visibility and map readability.
• Bing Map Access: If type = "bing" fails, you may need to set a Bing Maps API key with options(openstreetmap.bing.key = "your_key"), or switch to type = "osm" for OpenStreetMap instead.

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