# Find if point is within geojson polygons

**URL:** <https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416>\
**Category:** Geo\
**Tags:** polygons, archgdal, libgeos\
**Created:** [October 26, 2021, 4:43pm UTC](https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416 "2021-10-26T16:43:44Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [October 26, 2021, 4:43pm UTC](https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416/1 "2021-10-26T16:43:44Z")

</div>

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:

```julia

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

```julia
p = lat:70, lon:-40,

```

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

Any help would be greatly appreciated.

---

<div class="post-metadata">

**Author:** ![evanfields](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evanfields/32/1744_2.png) [@evanfields](https://discourse.julialang.org/u/evanfields)\
**Post date:** [October 26, 2021, 7:19pm UTC](https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416/2 "2021-10-26T19:19:48Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [October 26, 2021, 9:31pm UTC](https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416/3 "2021-10-26T21:31:24Z")

</div>

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](https://github.com/openlayers/openlayers/blob/dc8c9c/examples/data/geojson/polygon-samples.geojson)

```julia
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:

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

```

---

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [October 26, 2021, 10:47pm UTC](https://discourse.julialang.org/t/find-if-point-is-within-geojson-polygons/70416/4 "2021-10-26T22:47:36Z")

</div>

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

```julia
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:

```julia
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!