Find if point is within geojson polygons

Specific Domains Geo
, ,
1

I am slowly dipping my toe into the Julia waters and I’m hoping someone could point me in the right direction as to how I can determine if a point falls within any of the polygons defined by a geojson file. I’ve found GeoJSON’s dict2geo for converting from a dictionary to a geointerface:


cubejson_path = "/its-live-data.jpl.nasa.gov/datacubes/v01/datacubes_100km_v01.json"
cubejson = AWS.AWSServices.s3("GET", cubejson_path)
cubes = dict2geo(cubejson)

but I’m struggling to define a single point

p = lat:70, lon:-40,

and determine if falls within a any of the feature polygons.

Any help would be greatly appreciated.

2

I’m not a geospatial expert, so anyone should correct errors in the following, but:
I believe GeoJSON.jl and GeoInterface.jl don’t provide direct geometric operations on geometry objects. To test if a lat-lon point is within a polygon, you’ll probably want to use LibGEOS.jl, which wraps the GEOS library for planar geometry.

LibGEOS objects can be instantiated from similarly typed objects (eg Point -> LibGEOS.Point) GeoInterface objects, and LibGEOS offers the geometric predicates you might want like overlaps.

1 Like
3

I cannot access the linked data, probably I need to setup access tokens. But here is an example of how you could do this with ArchGDAL. With it you can directly open data from S3 (through curl), and do point in polygon operations (through GEOS).

This is a polygon GeoJSON file we can use as an example: https://github.com/openlayers/openlayers/blob/dc8c9c/examples/data/geojson/polygon-samples.geojson

using ArchGDAL, DataFrames

# prepend /vsicurl/ to read directly through curl, see https://gdal.org/user/virtual_file_systems.html
path = "/vsicurl/https://raw.githubusercontent.com/openlayers/openlayers/dc8c9cfabb47790fe404cf9156d281e0a8f38ae3/examples/data/geojson/polygon-samples.geojson"
dataset = ArchGDAL.read(path)
layer = ArchGDAL.getlayer(dataset, 0)
df = DataFrame(layer)

point = ArchGDAL.createpoint(-71.54, 47.60)

for row in eachrow(df)
    if ArchGDAL.contains(row[1], point)
        println("point in $(row.name)")
    end
end
# => point in L'Étoile-du-Nord

I’m guessing if your authentication is set up correctly this must work:

cubejson_path = "/vsis3/its-live-data.jpl.nasa.gov/datacubes/v01/datacubes_100km_v01.json"
dataset = ArchGDAL.read(cubejson_path)
2 Likes
4

@visr and @evanfields thank you both for the guidance. @visr had it exactly correct. I couldn’t get public s3 read to work:

AWS.global_aws_config(AWSConfig(creds=nothing, region = "us-west-2"))
path2json = "s3://its-live-data.jpl.nasa.gov/datacubes/v01/datacubes_100km_v01.json"
dataset = ArchGDAL.read(path2json)

But I was successful using the url:

path2json = "http://its-live-data.jpl.nasa.gov.s3.amazonaws.com/datacubes/v01/datacubes_100km_v01.json"
dataset = ArchGDAL.read(path2json)
layer = ArchGDAL.getlayer(dataset, 0)
df = DataFrame(layer)

point = ArchGDAL.createpoint(47.60, -71.54)

for row in eachrow(df)
    if ArchGDAL.contains(row[1], point)
        println("point in $(row.data_epsg)")
    end
end

U both rock, thanks!

1 Like