Speed up ArchGDAL.contains for thousands of data

Specific Domains Geo
1

Hi Guys,
I have geopoints like this

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

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

and I got a function

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!

2

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

3

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

#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

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!

4

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

out =  gmtread(url);

or

out = gmtread("/vsicurl/" * url);
5

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

image
image1915×249 56.7 KB

6

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

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

7

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

8

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

9

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