特定区域卫星地图下载及R中SpatialPointsDataFrame分析相关技术咨询

Hey there! Great job getting that trjct_pts SpatialPointsDataFrame with your intensity sums sorted out. Let's walk through practical steps to build on this, whether you're looking to download the targeted satellite map, visualize your intensity data, or export your work.

First, let's recap your analysis output for clarity:

trjct_pts
class : SpatialPointsDataFrame
features : 104559
extent : 696621.4, 696688.7, 167659.2, 167739.8 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs
variables : 1
names : Intensity_kJ
min values : 161.951
max values : 1192.526

1. Download the Targeted Satellite Map Matching Your Extent #

Your data uses the Swiss Oblique Mercator projection (somerc), so first we'll convert it to WGS84 (EPSG:4326) since most satellite map services rely on this. We'll use the rosm and sp packages for this workflow:

# Install required packages if you haven't already
# install.packages(c("rosm", "sp", "ggplot2"))

library(sp)
library(rosm)
library(ggplot2)

# Convert your SpatialPoints to WGS84 projection
trjct_pts_wgs84 <- spTransform(trjct_pts, CRS("+init=epsg:4326"))

# Extract the geographic extent of your data to target the satellite map
map_extent <- bbox(trjct_pts_wgs84)

# Download satellite imagery tiles (adjust zoom for detail)
sat_map <- osm.raster(
  bbox = map_extent,
  type = "satellite",
  zoom = 18  # Higher zoom = more detail, adjust based on your area size
)

2. Visualize Your Intensity_kJ Data on the Satellite Map #

Combine your intensity points with the satellite base map to spot spatial patterns in your data. Here are two options:

Static Plot with ggplot2 #

# Convert SpatialPoints to a data frame for ggplot compatibility
trjct_df <- as.data.frame(trjct_pts_wgs84)
colnames(trjct_df)[1:2] <- c("lon", "lat")  # Rename coordinates for clarity

# Plot satellite base + intensity points
ggplot() +
  layer_spatial(sat_map) +  # Add satellite background
  geom_point(
    data = trjct_df,
    aes(x = lon, y = lat, color = Intensity_kJ),
    alpha = 0.6, size = 0.8
  ) +
  scale_color_viridis_c(option = "plasma", name = "Intensity (kJ)") +
  theme_minimal() +
  labs(title = "Intensity Distribution Over Satellite Map")

Interactive Map with leaflet #

library(leaflet)

leaflet(trjct_pts_wgs84) %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
  addCircleMarkers(
    color = ~colorNumeric("plasma", Intensity_kJ)(Intensity_kJ),
    radius = 2, opacity = 0.7,
    popup = ~paste("Intensity:", round(Intensity_kJ, 2), "kJ")
  ) %>%
  addLegend(pal = colorNumeric("plasma", trjct_pts_wgs84$Intensity_kJ), 
            values = ~Intensity_kJ, title = "Intensity (kJ)") %>%
  addLayersControl(baseGroups = c("Satellite"))

3. Export Your Results #

• Export the SpatialPointsDataFrame to a shapefile (for use in GIS tools like QGIS):

  library(rgdal)
  writeOGR(trjct_pts, dsn = "./intensity_points", layer = "trjct_intensity", driver = "ESRI Shapefile")
  

• Export your static map as a high-resolution image:

  ggsave("intensity_satellite_map.png", width = 10, height = 8, dpi = 300)
  

If you're aiming for something specific—like rasterizing your intensity data into a heatmap, or customizing the satellite map styling—just let me know and I can refine these examples further!

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