# Speed up ArchGDAL.contains for thousands of data

**URL:** <https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847>\
**Category:** Geo\
**Created:** [March 14, 2022, 6:29am UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847 "2022-03-14T06:29:45Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 14, 2022, 6:29am UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/1 "2022-03-14T06:29:45Z")

</div>

Hi Guys,  
I have geopoints like this

```nohighlight
import ArchGDAL as AG
lon = LinRange(-180,180,7500)
lat = LinRange(-90,90,5500)

function meshgrid(x, y)
    X = [i for i in x, j in 1:length(y)]
    Y = [j for i in 1:length(x), j in y]
    return X, Y
end

Lon, Lat = meshgrid(lon,lat)
geolist = AG.createpoint.(reshape(Lon,1,:),reshape(Lat,1,:))
geolist

```

and a DataFrame

```nohighlight
data = DataFrame(AG.getlayer(AG.read("xxxx.geojson"),0))

```

and I got a function

```nohighlight
function index2score(index)
    temp_score = sum([AG.contains(tmp[index,1], xloc) for xloc in geolist])
    #I don't know why this one could not work
    #temp_score = sum(AG.contains.(tmp[index,1],geolist))
    name = tmp[index, "name"]
    println("name is $name and score is $temp_score")
    #return DataFrame(name = [name], score = [tmp_score])
end

```

It took me 20 minutes to run this, is there any way to speed up?  
PS: I am a newbie to use ArchGDAL, I get I could get more efficiency there  
Thanks!

---

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 14, 2022, 7:43am UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/2 "2022-03-14T07:43:16Z")

</div>

Sorry guys, it seems to be a replicated question, I would Try using `GMT.jl` first.

> [@Efficiently check if points are contained in polygons](https://discourse.julialang.org/t/efficiently-check-if-points-are-contained-in-polygons/74415/5):
>
> Thank you for your reply! I will change the array data type. Thank you for pointing out these C packages, do you know where to find tutorials for them? For now I will focus on learning Julia, but it would be nice to have a look at these packages.

---

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 14, 2022, 11:47am UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/3 "2022-03-14T11:47:47Z")

</div>

It seems I have to convert the `DataFrame` to `geojson` type to use `gmtread` function  
Here is my function to get the data

```nohighlight
#This one use ArchGDAL and it works
function get_geodata(url::String)
    final_body = String(HTTP.get(url).body) # this should be a geojson file
    geo_file = AG.read(final_body) 
    geo_data = DataFrame(AG.getlayer(geo_file,0))
    return geo_data
end
## This one does not work (the kernel would die out immediately)
function gmt_geodata(url::String)
    final_body = String(HTTP.get(url).body)
    return gmtread(final_body;dataset=true)
end

```

and the `download_country` function would require to call the `get_geodata` function several times and then `vcat` all `DataFrame`

```nohighlight
function download_country(cnmap::chinamap,target="边界")
        condition = (cnmap.raw_data[:,"adcode_third"] .== "00") .&& (cnmap.raw_data[:,"adcode_second"] .== "00") .&& !(cnmap.raw_data[:,"name"] in (["x1", "x2"]))
        tmp_data = cnmap.raw_data[condition,1:end]
        geo_province_datas = [download_province(cnmap,i,target) for i in tmp_data.name] # show return a vector of DataFrame
        return vcat(geo_province_datas...)

```

So any suggestions to make GMT work in this case!! Really appreciate all you guys!

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 14, 2022, 12:11pm UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/4 "2022-03-14T12:11:34Z")

</div>

I don’t have time to look at this before the end of the day but try just this

```julia
out = gmtread(url);

```

or

```julia
out = gmtread("/vsicurl/" * url);

```

---

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 14, 2022, 12:26pm UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/5 "2022-03-14T12:26:36Z")

</div>

Thanks, but it seems not working.  
How come the return to be a vector of GMTdataset type

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/a/7a7f5359880a2f1e33a2230b2908c341b7d7c3de.png)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 14, 2022, 6:48pm UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/6 "2022-03-14T18:48:23Z")

</div>

> [@Frankiewaang](#):
>
> Thanks, but it seems not working.

Why not? It failed when apparently gave ii a WMs or MFS url, but worked when passing in a json file.

> [@Frankiewaang](#):
>
> How come the return to be a vector of GMTdataset type

Because that’s what it should return when reading vector data. What were you expecting?

---

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 15, 2022, 3:35am UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/7 "2022-03-15T03:35:25Z")

</div>

I get it. Thanks for your explanation. I would look for another way to solve it.

---

<div class="post-metadata">

**Author:** ![ahmoreira](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmoreira/32/32661_2.png) [@ahmoreira](https://discourse.julialang.org/u/ahmoreira)\
**Post date:** [March 18, 2022, 4:01pm UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/8 "2022-03-18T16:01:08Z")

</div>

@Frankiewaang ,just to let you know that after testing these functions with a larger dataset (40.000 points) inpolygon from Rasters.jl performed better.

---

<div class="post-metadata">

**Author:** ![Frankiewaang](https://avatars.discourse-cdn.com/v4/letter/f/f9ae1b/32.png) [@Frankiewaang](https://discourse.julialang.org/u/Frankiewaang)\
**Post date:** [March 25, 2022, 12:39pm UTC](https://discourse.julialang.org/t/speed-up-archgdal-contains-for-thousands-of-data/77847/9 "2022-03-25T12:39:52Z")

</div>

Thanks! Your thread helps me a lot. But I can’t get `Rasters.jl` to read `geojson` file, and it seems there are no solutions but to convert it to `shp` file