如何判断经纬度是否属于多边形并处理zone文件中的坐标字符串

Got it, let's break down how to solve this problem—you need to check if a given lat-lon point falls into any of the 1994 polygons from your zone data file. Since your sample syntax looks like R, I'll walk through a complete R-based solution step by step:

First, we need to read and clean up the raw data into a usable format. Each line in your file has a zone ID and a string of polygon vertices. Here's how to parse it:

# Read the zone data file (adjust the file path as needed)
zone_lines <- readLines("zone.txt")

# Parse each line into zone ID and structured polygon vertices
zone_polygons <- lapply(zone_lines, function(line) {
  # Split the line into ID and coordinate parts (matches your sample format)
  parts <- strsplit(line, '", "')[[1]]
  zone_id <- gsub('^"', '', parts[1])  # Remove leading quote from ID
  coord_str <- gsub('"$', '', parts[2]) # Remove trailing quote from coordinates
  
  # Split coordinate string into individual vertex pairs
  vertex_strings <- strsplit(coord_str, ', ')[[1]]
  
  # Convert each vertex to numeric lat-lon values
  vertices <- do.call(rbind, lapply(vertex_strings, function(v) {
    as.numeric(strsplit(v, ' ')[[1]])
  }))
  colnames(vertices) <- c("lat", "lon")
  
  # Return a list with zone ID and polygon vertices
  list(zone_id = zone_id, polygon = vertices)
})

The ray casting algorithm is a reliable way to check if a point is inside a polygon. Here's a custom implementation:

# Ray casting function to determine if a point is inside a polygon
point_in_polygon <- function(target_point, polygon) {
  # target_point: vector in format c(lat, lon)
  # polygon: n x 2 matrix of lat-lon vertices
  
  n <- nrow(polygon)
  # Close the polygon if the first and last vertices don't match
  if (!all(polygon[1, ] == polygon[n, ])) {
    polygon <- rbind(polygon, polygon[1, ])
    n <- n + 1
  }
  
  inside <- FALSE
  for (i in 1:(n - 1)) {
    p1 <- polygon[i, ]
    p2 <- polygon[i + 1, ]
    
    # Check if the point's latitude falls within the edge's latitude range
    lat_in_range <- ((p1[1] > target_point[1]) != (p2[1] > target_point[1]))
    # Check if the point's longitude is left of the edge's intersection with the point's latitude
    lon_intersect <- (target_point[2] < (p2[2] - p1[2]) * (target_point[1] - p1[1]) / 
                        (p2[1] - p1[1]) + p1[2])
    
    if (lat_in_range && lon_intersect) {
      inside <- !inside
    }
  }
  inside
}

Use the Function to Check a Point #

Let's test with an example point:

# Example target point (lat=25, lon=35)
target_point <- c(25, 35)

# Find all zones that contain the point
matching_zones <- lapply(zone_polygons, function(zone) {
  if (point_in_polygon(target_point, zone$polygon)) {
    zone$zone_id
  } else {
    NULL
  }
})
# Remove empty entries and print results
matching_zones <- unlist(matching_zones)
cat("The point belongs to zone(s):", paste(matching_zones, collapse = ", "), "\n")

For large datasets (1994 polygons), using a dedicated spatial package like sf is more efficient and less error-prone. Here's how to use it:

# Install and load sf if not already installed
if (!require(sf)) {
  install.packages("sf")
  library(sf)
}

# Convert parsed polygons to sf spatial objects
# Note: sf uses lon-lat order, so we swap the lat/lon columns
sf_polygons <- lapply(zone_polygons, function(zone) {
  # Create a polygon geometry (sf expects vertices in lon-lat order)
  polygon_geom <- st_polygon(list(matrix(c(zone$polygon[, "lon"], zone$polygon[, "lat"]), ncol = 2)))
  # Create an sf object with zone ID and geometry (WGS84 coordinate system)
  st_sf(zone_id = zone$zone_id, geometry = polygon_geom, crs = 4326)
})
# Combine all polygons into a single sf object
sf_zones <- do.call(rbind, sf_polygons)

# Convert the target point to an sf point (again, lon-lat order)
target_sf <- st_sfc(st_point(c(target_point[2], target_point[1])), crs = 4326)

# Perform the point-in-polygon intersection check
matches <- st_intersects(target_sf, sf_zones, sparse = FALSE)
matching_zones_sf <- sf_zones$zone_id[matches[1, ]]

cat("Using sf, the point belongs to zone(s):", paste(matching_zones_sf, collapse = ", "), "\n")

Key Notes #

• Coordinate Order: Most spatial libraries (including sf) use lon-lat order, but your data uses lat-lon. Make sure to swap columns when working with sf to avoid incorrect results.
• Polygon Closure: The ray casting function automatically closes polygons if the first and last vertices don't match—this is required for the algorithm to work correctly.
• Data Consistency: Ensure your zone file has consistent formatting (each line follows the "ID", "coord1, coord2,..." structure). If your actual file differs, adjust the strsplit parameters in the parsing step.

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