Amy Xiejing
Amy Xiejing

Reputation: 23

Filter polygons on coordinates

I'm currently working with a shapefile that does not contain explicit latitude and longitude coordinates. My goal is to filter out areas with latitudes above 60 degrees. I want to delete the small area in the lower right corner, but I don't have the latitude and longitude. Can anyone suggest a code to achieve this?

The shapefile look like this:

enter image description here

This is my current code:

shapefile_path <- "/Users/langbai/Desktop/lcsd000b16a_e/lcsd000b16a_e.shp"

shapefile_data <- st_read(shapefile_path)

convex_hulls_inset <- convex_hulls_with_farms %>%
  filter(lat < 60)

shapefile_inset <- shapefile_data %>%
  filter(st_coordinates(st_centroid(geometry))[, 2] < 80)

Upvotes: 1

Views: 202

Answers (2)

Spacedman
Spacedman

Reputation: 94182

Here's how to select based on any point in the polygon being above the threshold latitude.

Sample data:

library(sf)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))

Get full matrix of X, Y and L1, L2, L3 components:

xy = st_coordinates(nc)

Now select rows of that matrix where the Y column is above the threshold:

thresh=35
above = xy[ xy[,"Y"] > thresh, ]

Then the unique values of L3 are the rows that are on or above the threshold:

row.above = unique(above[,"L3"])

And we can plot them with a line at the threshold:

plot(nc$geometry)
plot(nc$geometry[row.above], col="red", add=TRUE)
plot(nc$geometry[-row.above], col="blue", add=TRUE)
abline(h=thresh, lwd=3)

enter image description here

You can adjust the logic depending on if you want features that are completely above or completely below the threshold:

below = xy [ xy[,"Y"] < thresh, ]
row.below = unique(below[,"L3"])
plot(nc$geometry)
plot(nc$geometry[row.below], col="blue", add=TRUE)
plot(nc$geometry[-row.below], col="red", add=TRUE)
abline(h=thresh, lwd=3)

enter image description here

Upvotes: 0

L Tyrone
L Tyrone

Reputation: 6885

First, create a new column using mutate(), then filter() that column. Here's an example using the built-in nc dataset:

library(sf)
library(dplyr)
library(ggplot2)

# Example data
nc <- st_read(system.file("shape/nc.shp", package = "sf"))

shapefile_inset <- nc %>%
  mutate(y = st_coordinates(st_centroid(geometry))[, 2]) %>%
  filter(y < 35) %>%
  select(-y)

ggplot() + 
  geom_sf(data = nc) +
  geom_sf(data = shapefile_inset, fill = "#DF536B")

result

Upvotes: 1

Related Questions