I want to know if a lat/long coordinate is inside a specified set of US states. I downloaded a shapefile for the US states from the Census Bureau. Suppose my coordinate is (38.4511,-80.5194), and I want to know if this coordinate is in the states of PA, NY, or NJ.
My current code is the following:
using Shapefile, GeoInterface, GeometryBasics
state_table = Shapefile.Table("path/to/cb_2018_us_state_500k/cb_2018_us_state_500k.shp")
## data frame for selecting states
state_DF = DataFrame(state_table)
keep_states = ["NJ", "PA", "NY"]
state_DF = filter(:STUSPS => x -> x in keep_states, state_DF)
function is_in(point::GeometryBasics.Point, polygon)
return in(point, polygon)
end
for i in 1:3
if is_in(Point(-80.5194, 38.4511), state_DF[i,1])
println("The point is inside ", state.attributes["STATE_NAME"])
break # Stop after finding the first matching state
end
end
Then, I get the following error:
ERROR: MethodError: no method matching iterate(::Shapefile.Polygon)
Closest candidates are:
iterate(::Union{LinRange, StepRangeLen}) at range.jl:872
iterate(::Union{LinRange, StepRangeLen}, ::Integer) at range.jl:872
iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at dict.jl:712
...
How can I accomplish my task? I’ve tried bypassing this error for a bit now.