[QST] rasterization solution for data tiles using PIP?

[gouldingken]

After attending @thomcom's FOSS4G presentation, it occured to me that I have a use case for your point-in-polygon (PIP) algorithm that I would never have thought to explore through something like this!

I have a bit of a crazy use case - which is rasterizing the entire US at pretty high res for census blocks (slippy map zoom level 13 see zoom tile reference). This is the first step in distributing block-level census data onto those tiles in an easy to query format for geo-spatial viz.

The methodology as I imagine it:

• For each 256x256 tile, I'd load in the vector polygon geometries, then test point-in-polygon using your algorithm applied to a single test point for every pixel (actually 255x255 points because of how they fall on the grid).
• Each pixel would be assigned the ID of the geometry it's in (there are no overlapping geometries) or NULL if there is no containing polygon.

Some stats:

• I have 64k zoom 11 vector tiles (about 250MB of compressed PBF data)
• each vector tile will be rasterized into 16 256x256 tiles at zoom 13 (giving about a million raster tiles).
• the total number of PIP operations is around 67 billion (65k pixels per 256x256 tile)
• final output will be in geo-png-db format, but can be easily converted after-the-fact

I'd love to find a way to handle this on a powerful GPU! I've been trying to use parallelized serverless node scripts to run PIP, but it's pretty sluggish. I had initially tried graphical rendering methods for the rasterization, but ended up brute-forcing it with PIP because the results are better when polygons get small, and where geometries meet.

Census block polygons are vector tiles (currently geobuf pbf format). The tile pixel points can be created as a point-data geospatial file if easier - or potentially calculated on the fly.

I'd appreciate any advice you could offer regarding this use case,

Thanks,
Ken

[harrism]

GPUs have very fast hardware rasterizers, so it seems like a shame to implement rasterization using per-pixel point-in-polygon tests. That said, it would require using a graphics API to access hardware rasterization, which may not be something we want to do in cuSpatial. Still, my gut tells me there is a better way than PIP.

[thomcom]

Hey Ken, do you have an update on your progress here? I've had some ideas I'd like to chat about.

[gouldingken]

Hi Thomson, the project that requires this 'rasterization' process keeps getting bumped by more pressing tasks, so we haven't made any progress on our end. We're still really interested in how to do this however and I'd love to hear your ideas. What's the best forum for a discussion?

[gouldingken]

I didn't respond to harrism's suggestion. I know nothing of hardware rasterization, but I will say that all the software rasterization techniques we tried fell short in terms of accurately representing adjacent polygons and leaving no gaps. PIP is obviously much slower, but gives the results we need.

[harrism]

Interesting. Did you try cuSpatial's spatial join for this? There is room for optimization since your data are so regular (pixels) but you could generate points from the pixels (not all 67 billion at once, but you could fit O(1e9) on a GPU at a time) and index them into a quadtree and then query the polygon geometry against all 1e9 points at once. On an A100 GPU we achieve over 100B PIP queries per second including the time to build the quadtree (this is for a dataset with 168M points and 263 polygons -- performance is somewhat dependent on polygon complexity and number of polygons). Extrapolating, doing this job in less than a second looks possible with existing APIs.

Due to the regularity of the points (regular grid), the quadtree will be very regular as well. So it might be worth exploring shortcuts in the PIP algorithm for this use case (e.g. faster building of the quadtree since it is fully regular, and then rasterizing entire contained quadrants).