Mencari solusi tercepat untuk analisis Point in Polygon 200 juta poin [ditutup]

Saya memiliki CSV yang berisi 200 juta pengamatan dengan format berikut:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"

Untuk setiap set koordinat (x1 / y1 dan x2 / y2), saya ingin menetapkan Blok Sensus AS atau Blok Sensus yang termasuk di dalamnya (Saya mengunduh traktus TIGER traktus sensus TIGER di sini: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Jadi, saya perlu melakukan operasi point-in-poligon dua kali untuk setiap pengamatan. Sangat penting bahwa pertandingan sangat akurat.

Apa cara tercepat untuk melakukan ini, termasuk waktu untuk mempelajari perangkat lunak? Saya memiliki akses ke komputer dengan Memori 48GB - dalam kasus yang mungkin merupakan kendala yang relevan.

Beberapa utas merekomendasikan penggunaan PostGIS atau Spatialite (Spatialite terlihat lebih mudah digunakan - tetapi apakah ini seefisien PostGIS?). Jika itu adalah pilihan terbaik, apakah penting untuk mengisi Indeks Spasial (RTree?)? Jika demikian, bagaimana cara melakukannya (mis. Menggunakan Census Tract Shapefile)? Saya akan sangat berterima kasih atas rekomendasi yang menyertakan kode contoh (atau pointer ke kode contoh).

Upaya pertama saya (sebelum menemukan situs ini) terdiri dari menggunakan ArcGIS untuk melakukan join spasial (hanya x1 / y1) dari subsampel data (100.000 poin) pada Blok Sensus AS. Butuh waktu lebih dari 5 jam sebelum saya membunuh prosesnya. Saya berharap untuk solusi yang dapat diterapkan pada seluruh dataset dalam waktu komputasi kurang dari 40 jam.

Permintaan maaf karena mengajukan pertanyaan yang telah diajukan sebelumnya - Saya sudah membaca jawaban, dan saya bertanya-tanya bagaimana menerapkan rekomendasi. Saya tidak pernah menggunakan SQL, Python, C, dan hanya pernah menggunakan ArcGIS sebelumnya - saya seorang pemula yang lengkap.

ST_DWithin lebih cepat dalam pengujian saya daripada ST_Intersects. Itu mengejutkan, terutama karena algoritma geometri yang sudah disiapkan seharusnya digunakan untuk kasus-kasus seperti ini. Saya pikir ada kemungkinan ini akan jauh lebih cepat daripada yang saya tunjukkan di sini.

Saya melakukan beberapa tes lagi dan dua hal kecepatannya hampir dua kali lipat. Pertama, saya mencoba di komputer yang lebih baru, tetapi masih laptop yang cukup biasa, mungkin kecuali dari SATA3 ssd -disks.

Kemudian permintaan di bawah ini memakan waktu 18 detik, bukan 62 detik di laptop lama. Selanjutnya saya menemukan bahwa saya benar-benar salah sebelumnya ketika saya menulis bahwa indeks pada tabel-titik tidak perlu. Dengan indeks itu, ST_Intersects bertindak sesuai yang diharapkan dan segalanya menjadi sangat cepat. Saya meningkatkan jumlah poin di tabel poin menjadi 1 juta poin dan permintaan:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct , t WHERE ST_Intersects(imported_ct.geom , t.geom);

berjalan dalam 72 detik. Karena ada 1249 poligon, uji 1249000000 dilakukan dalam 72 detik. Itu menghasilkan sekitar 17000000 tes per detik. Atau menguji hampir 14.000 poin terhadap semua poligon per detik.

Dari pengujian ini, 400000000 poin Anda untuk diuji harus memakan waktu sekitar 8 jam tanpa kesulitan dengan mendistribusikan beban ke beberapa core. PostGIS tidak pernah berhenti mengesankan saya :-)

Pertama, untuk memvisualisasikan hasil, Anda dapat menambahkan titik geometri ke tabel yang dihasilkan, buka di QGIS misalnya dan gaya dengan nilai unik pada bidang import_ct.

Kedua, ya, Anda juga bisa mendapatkan poin yang berada di luar poligon apa pun dengan menggunakan gabungan kanan (atau kiri) seperti ini:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct right join t ON ST_Intersects(imported_ct.the_geom , t.geom);

Saya melakukan beberapa tes untuk memverifikasi jika tampaknya mungkin PostGIS.

Pertama, sesuatu yang tidak saya mengerti. Anda memiliki dua poin per baris. Apakah kedua poin selalu dalam poligon yang sama? Maka cukup melakukan perhitungan pada salah satu poin. Jika mereka bisa dalam dua poligon yang berbeda, Anda akan membutuhkan cara untuk menghubungkan satu baris baris ke dua poligon.

Dari tes ini sepertinya bisa dilakukan, tetapi Anda mungkin memerlukan beberapa solusi kreatif untuk menyebarkan beban lebih dari satu cpu-core.

Saya menguji pada laptop berusia 4 tahun dengan cpu dual core centrino (kira-kira 2.2GHz saya kira), RAM 2GB. Jika Anda memiliki 48 BG RAM, saya kira Anda juga memiliki lebih banyak cpu-power.

Apa yang saya lakukan adalah membuat tabel poin acak dengan 100000 poin seperti ini:

CREATE TABLE t AS
WITH r AS
(SELECT ST_Extent(the_geom)::geometry ext FROM imported_ct)
SELECT ST_Point(x,y) AS geom FROM 
(SELECT GENERATE_SERIES(1,100000)) s,
(SELECT ST_Xmin(ext)+(random()*(ST_Xmax(ext)-ST_Xmin(ext))) x, ST_Ymin(ext)+(random()*(ST_Ymax(ext)-ST_Ymin(ext))) y FROM r
) f;

Kemudian menambahkan gid seperti:

ALTER TABLE t ADD COLUMN GID SERIAL;

Kemudian jalankan:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE ST_Dwithin(imported_ct.the_geom , t.geom,0);

membutuhkan sekitar 62 detik (bandingkan dengan hasil ArcGIS Anda dengan jumlah poin yang sama). Hasilnya adalah tabel yang menghubungkan titik-titik di tabel saya t dengan gid di tabel dengan saluran sensus.

Dengan kecepatan itu Anda akan melakukan 200 poin pabrik dalam waktu sekitar 34 jam. Jadi, jika cukup dengan memeriksa salah satu intinya, laptop lama saya dapat melakukannya dengan satu inti.

Tetapi jika Anda perlu memeriksa kedua poin itu mungkin lebih sulit.

Kemudian Anda dapat secara manual mendistribusikan beban ke lebih dari satu inti dengan memulai beberapa sesi melawan db dan menjalankan kueri yang berbeda.

Dalam contoh saya dengan 50.000 poin dan dua cpu-core saya mencoba:

CREATE TABLE t1 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid >50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

pada satu db-sesi sekaligus dengan menjalankan:

CREATE TABLE t2 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid <=50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

pada sesi db lain.

Itu memakan waktu sekitar 36 detik sehingga sedikit lebih lambat dari contoh pertama mungkin tergantung pada penulisan disk pada saat yang sama. Tetapi karena core bith bekerja pada saat yang sama, tidak memerlukan lebih dari 36 detik waktu saya.

Untuk tabel gabungan t1 dan t2 a dicoba:

CREATE TABLE t3 AS 
SELECT * FROM t1
UNION ALL
SELECT * FROM t2;

menggunakan sekitar setengah detik.

Jadi, dengan perangkat keras yang lebih segar dan mendistribusikan beban ke banyak core, ini harus benar-benar mungkin bahkan jika dunia nyata akan lebih lambat daripada test case.

Mungkin perlu dicatat bahwa contohnya adalah dari Linux (Ubuntu). Menggunakan Windows akan menjadi cerita lain. Tetapi saya memiliki semua aplikasi harian lainnya berjalan sehingga laptop ini cukup banyak dimuat dari sebelumnya. Jadi itu mungkin mensimulasikan case windows cukup baik mungkin, tanpa membuka apa pun kecuali pgadmin.

Mungkin cara termudah adalah dengan PostGIS. Ada beberapa tutorial di internet cara mengimpor data titik csv / txt ke PostGIS. Tautan1

Saya tidak yakin tentang kinerja pencarian titik-dalam-poligon di PostGIS; itu harus lebih cepat dari ArcGIS. Indeks spasial GIST yang digunakan PostGIS cukup cepat. Link2 Link3

Anda juga bisa menguji indeks geospasial MongoDB . Tetapi ini membutuhkan sedikit lebih banyak waktu untuk memulai. Saya percaya bahwa MongoDB bisa sangat cepat. Saya belum mengujinya dengan pencarian titik-dalam-poligon jadi tidak bisa memastikan.