3

I have no programming skills, but am thinking I may need to at least figure out how to run this one script from inside QGIS.

I have a vector file of global quarter degree grid cells and a point file containing an incredibly large number of points with population values.

I want to create a polygon file that has the sum of the population values (not just the count of the points) for each grid cell.

This is really just a basic Join by Location (Spatial Join), but the files are so large it will not work.

Can anyone help me?

  • What data format are the layers (shapefile, sqlite, PostGIS...)? – artwork21 Oct 4 '17 at 16:57
  • I believe they are shapefiles. One is a vector grid I made in QGIS and the other is a point file I generated in ArcMap – Arthur_15 Oct 4 '17 at 17:07
  • I love that people are viewing this and commenting. gives me hope that someone might have an answer for me :) – Arthur_15 Oct 4 '17 at 19:46
  • Maybe convert it into a sqlite/spatialite db or PostgreSQL/PostGIS layer where you can take advantage of spatial indexes or split out your heavy data into smaller chunks/layers. – artwork21 Oct 4 '17 at 19:53
2

I find for your particular scenario - points versus a regular grid - I go for a programatic approach (e.g. Python).

I do this for a couple of reasons:-

  • it's faster than a traditional spatial join, as it's much easier to test for intersections with cells.
  • unlike a database lookup or spatial join, it's easier to monitor progress :)

Also, there's no need to consider every 1/4 degree cell if the data is 'sparse'. The world has 1,036,800 of those cells. 70% + of the cells will be empty (sea!) and most of the remaining land cells will be empty, too. That's a lot of wasted effort, however you do this - whether you use shapefiles or a database.

This code uses a CSV file for the points. I used the NASA Meteorite Database (from here) if you want to try this code out. In my case I'm summing up the weight of meteorites which land in each quarter-degree grid cell.

To keep things simple, it's standalone python and doesn't need any special libraries.

import csv

sums = {} # key=(cell_x,cell_y), value=sum

# load points file, this is a CSV (comma-separated) with long, lat and values in each record

with open("/tmp/Meteorite_Landings.csv","r") as fi:
    reader = csv.reader(fi,delimiter=",")
    count = 0
    for line in reader:
        try:
            lat = float(line[7]) # 8th column in my case for latitude
            lon = float(line[8]) # 9th column in my case for longitude
            val = float(line[4]) # value to sum is in 5th column
            cell_x = int((lon + 180.0) * 4) # 4 cells per degree
            cell_y = int((lat + 90.0) * 4)
            key = (cell_x, cell_y)
            if not key in sums:
                sums[key] = val
            else:
                sums[key] += val
        except:
            pass # ignore header line, lines with errors

# write this out to a TAB separated file.
# can load this into QGIS as a delimited layer.

with open("/tmp/foo.csv","w") as fo:
    fo.write("GEOM\tSUM\n")
    for key in sums:
        cell_x, cell_y = key
        total = sums[key]
        # convert cell number back to lat,lon of SW corner of cell
        nw_x = sw_x = -180.0 + (0.25 * cell_x) # convert to longitude of SW corner
        se_y = sw_y = -90.0 + (0.25 * cell_y)  # convert to longitude of SW corner
        ne_x = se_x = sw_x + 0.25
        nw_y = ne_y = sw_y + 0.25
        poly = "POLYGON(({} {},{} {},{} {},{} {},{} {}))".format(sw_x,sw_y,nw_x,nw_y,ne_x,ne_y,se_x,se_y,sw_x,sw_y)
        fo.write("{}\t{}\n".format(poly, total))

print("Finished!")

The output is CSV which you can then import into QGIS as a delimited layer. Only those cells with a value exist using this approach..

enter image description here

This summed up 30k records against a 1/4 degree grid in under a second, on a 10-year old laptop.

CSV isn't exactly the best format for huge amounts of data, but I'm guessing your data is sparse enough that this doesn't matter.

1

I'd try a QGIS 3.0 nightly release first - Select/Extract by Location are many orders of magnitude faster in 3.0.

  • thanks for these answers. Excited to try it out. I'm skeptical of my ability to use Python, but maybe I will be able to figure it out. – Arthur_15 Oct 5 '17 at 16:05
  • do you have a link where I can download 3.0? – Arthur_15 Oct 5 '17 at 16:11
  • Is it actually 2.99 and is there an OS version? – Arthur_15 Oct 5 '17 at 16:29
  • Does anyone know where I can download QGIS 3.0? – Arthur_15 Oct 6 '17 at 19:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.