使用ppp函数绘制德州休斯顿炼油厂位置时出现绘图错误

Hey there! Let's work through that ppp() error you're facing while mapping refineries around Houston. First, let's break down the common pitfalls with spatstat's ppp() function and fix your workflow step by step.

Key Issues to Address #

ppp() requires two critical inputs: point data (your refinery locations) and a valid window (your study area boundary) formatted for spatstat. Your current code only sets up the county boundaries as SpatialPolygons objects, which aren't directly compatible with ppp(). Also, geographic coordinates (lat/lon) can cause errors since spatstat works best with planar projections (like UTM).

Step 1: Complete & Clean Your County Boundary Code #

First, let's finish the truncated chambers subset (Chambers County's GEOID is 48015) and merge your target counties:

# Load required packages
library(sp)
library(spatstat)
library(maptools)

# Define geographic projection (WGS84)
P4S.latlon <- CRS("+proj=longlat +datum=WGS84 ")
# Define planar projection (UTM Zone 10N, fits Houston area)
utm_proj <- CRS("+proj=utm +zone=10 +datum=WGS84 +units=m")

# Read county shapefile
county.lines <- readShapePoly("tl_2017_us_county/tl_2017_us_county.shp", verbose = T, proj4string = P4S.latlon)

# Extract target Houston-area counties
harris <- subset(county.lines, county.lines$GEOID == "48201")
brazoria <- subset(county.lines, county.lines$GEOID == "48039")
galveston <- subset(county.lines, county.lines$GEOID == "48167")
chambers <- subset(county.lines, county.lines$GEOID == "48015") # Complete this line

# Merge counties into a single study area
houston_counties <- rbind(harris, brazoria, galveston, chambers)

Step 2: Convert Boundaries to a spatstat Window #

ppp() won't accept raw SpatialPolygons—we need to convert them to an owin (observation window) object, and switch to a planar projection first:

# Convert county boundaries to UTM (planar projection)
houston_counties_utm <- spTransform(houston_counties, utm_proj)

# Convert to spatstat-compatible window
houston_window <- as.owin(houston_counties_utm)

Step 3: Prepare Refinery Points & Create ppp Object #

Assuming you have a data frame of refinery lat/lon coordinates, convert them to match the window's projection and build the ppp object:

# Example refinery data (replace with your actual data)
refineries <- data.frame(
  lon = c(-95.3698, -95.0122, -94.9036),
  lat = c(29.7604, 29.6456, 29.3847)
)

# Convert points to SpatialPoints and reproject to UTM
ref_sp <- SpatialPointsDataFrame(
  coords = cbind(refineries$lon, refineries$lat),
  data = refineries,
  proj4string = P4S.latlon
)
ref_sp_utm <- spTransform(ref_sp, utm_proj)

# Create the ppp object
ref_ppp <- ppp(
  x = coordinates(ref_sp_utm)[,1],
  y = coordinates(ref_sp_utm)[,2],
  window = houston_window
)

Step 4: Plot Your Data #

Now you can safely plot without errors:

plot(ref_ppp, main = "Refineries in Houston Metro Area", pch = 16, col = "darkred")
plot(houston_window, add = TRUE, border = "blue", lwd = 2)

Common Error Fixes #

If you're still hitting issues, check these:

• "Window must be an owin object": You forgot to convert your county polygons with as.owin().
• "Data are not in a plane": You're using lat/lon coordinates—always switch to a planar projection like UTM for spatstat.
• "Points outside window": Use clip.ppp(ref_ppp, houston_window) to remove points outside your study area.
• Deprecated readShapePoly: For a more modern workflow, use the sf package instead:

  library(sf)
  county.lines <- st_read("tl_2017_us_county/tl_2017_us_county.shp")
  houston_counties <- subset(county.lines, GEOID %in% c("48201", "48039", "48167", "48015"))
  houston_window <- as.owin(st_transform(houston_counties, utm_proj))
  

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