GRASS GSoC 2026 Parallelizing r.proj and Raster Processing Modules in GRASS

From GRASS-Wiki
Jump to navigation Jump to search
Student Name Kaushik Raja
Organization NumFOCUS
Mentor Name Huidae Cho, Anna Petrasova, Vaclav Petras
GitHub Fork View Repo
LinkedIn Profile View LinkedIn

Abstract

R.proj, r.param.scale, and r.geomorphon are some of the most commonly used modules in GRASS. R.proj reprojects raster maps between coordinate systems. This is important because real world data comes in many different projections, and maps have to be in the same one before they can be analyzed together. So many GRASS workflows that combine data sources start with r.proj. R.param.scale calculates terrain parameters like slope and curvature by sliding a window over an elevation map. R.geomorphon classifies every cell of an elevation map into a landform like a ridge, valley, or peak. It does this by looking outward from each cell in eight different directions and checks whether the terrain rises above or drops below the line of sight. That pattern of visible horizons tells it what shape the land around the cell is.

The problem with all three modules is that they have always been single threaded. On modern hardware most CPU cores sit idle while one core does all the work, so large maps take much longer time to work and return an output. The goal of this project is to parallelize them with OpenMP so that many threads work on the map at the same time.

The main challenge with parallelizing these modules is memory. To help manage memory and get the most speedup ratios possible, each module will work on a chunk of rows at a time, called a band. The band is sized so that the input it needs stays under a memory limit the user controls. The threads then split the rows of the band among themselves and work on them at the same time. In r.proj, each thread also gets its own file descriptor and its own PROJ object (for thread safety). Some projections bend so much that even one full row’s width is too much input to fit under the memory limit, so those are read in smaller column pieces instead. In r.param.scale the old sliding window was replaced with this same band design, and in between I found and fixed a data race in the GRASS math library.

As a result of this, users can use one of the three modules and get their outputs faster than before. Also, the band pattern used in these modules gives future contributors a template for parallelizing other raster modules.

Project Scope

  1. Parallelize r.param.scale with OpenMP under a user controlled memory limit
  2. Parallelize r.proj
  3. Parallelize r.geomorphon
  4. Add nprocs option so user can say how many threads form their machine they want to use
  5. Find and fix data races in the GRASS libraries that block modules from being parallelized
  6. Replace the old testsuites for r.param.scale, r.proj, and r.geomorphon with pytest tests
  7. Verify every parallel module produces output identical to the serial version
  8. Add a benchmarking script for each module and benchmark all modules across thread counts

Timeline

Period Timeline Tasks Status
Community Bonding Period May 1 - May 25
  1. Thread safety audit of gprojects library calls
  2. Study how the code in readcell.c works
  3. Set up benchmarking infrastructure and figure out how that works
  4. Look at remaining modules and understand how they work
  5. Finalize dev environment
  6. Confirm benchmarks are reproducible
  7. Agree on implementation details with mentors
Done
Official Coding Period May 25 - June 8
  1. Build r.proj proof of concept with RAM buffer and per thread PROJ contexts (PR #7185)
  2. Benchmark RAM buffer approach against tile cache approach
  3. Add user controlled memory option
  4. Build first r.param.scale parallel draft (PR #7236)
  5. Rule out a suspected Mac OpenMP bug, traced to benchmark setup
Done
June 9 - June 22
  1. Redesign r.param.scale from a single strip into a two level band and chunk structure that honors the memory option, modeled on r.neighbors (PR #7440)
  2. Give each thread its own input access so threads no longer share one reader
  3. Measure 4.6x speedup at 8 threads, up from 2.7x before the redesign, with output identical to the serial module
  4. Trace a scaling dip to thread load imbalance using per thread timing, correlation of 0.99 between imbalance and slow runs, ruling out a CPU frequency explanation
  5. Find a data race in G_ludcmp in the GRASS math library affecting every parallel caller, prove it with ThreadSanitizer and a controlled toggle experiment, report it (issue #7539)
  6. Fix the race and verify with ThreadSanitizer showing 14 race reports before and 0 after, repeated test runs going from 10 failures in 100 to 0 in 20, and a 40 configuration bit identical output comparison
  7. Check every caller of the fixed function at runtime to confirm safety for dependent modules
  8. Debug CI failures across platforms including CMake PROJ linkage and unguarded OpenMP timer calls, resolve a git branch divergence, unskip and pass the pytest suites on Python 3.10 and 3.13, pass all 26 checks, mark PR #7440 ready for review
Done
June 23 - July 6
  1. Design and build the band based r.proj (PR #7627)
  2. Size each band by projecting its edges back to the input to find the rows it needs
  3. Cut peak memory from 763 MB to 130 MB with identical output
  4. Measure 2.9x total and 5.2x compute speedup at 8 threads on a 105 million cell map
  5. Answer two rounds of mentor review
Done
July 7 - July 11
  1. Submit midterm evaluation
  2. Fix two CI failures, a missing PROJ dependency in the CMake build and OpenMP timer calls breaking the no OpenMP build
  3. Parallelize input reading with per thread file descriptors, 1.9x faster reads, about 3.0x total
  4. Move per thread PROJ context handling into the gproj library as new API functions at maintainer request
  5. Measure which projections need column splitting with a standalone footprint tool
  6. Split bands into column sections so tilted reprojections that failed under the memory cap now complete, verified identical on three projections
  7. Cut the section sizing search from 54 seconds to 2.6 seconds
  8. Benchmark everything and present results to mentors
Done
July 20 - August 3
  1. Implement the two speedup items from review on PR #7627, keeping the input strip resident across bands and writing the previous band's output while the next band computes
  2. Run bilinear, bicubic and lanczos through the banded parallel path
  3. Move the method reference tests into their own PR #7766, replacing the old testsuite with pytest
  4. Fix a data race on two globals in the projection library (PR #7764, merged)
  5. Run the full benchmark grid, post results and scaling graphs, mark PR #7627 ready for review
  6. Look at r.geomorphon and fix an already existing crash on regions smaller than the search window (PR #7773, merged)
  7. Parallelize r.geomorphon with per band strips, per thread file descriptors and a nprocs option (PR #7783), about 5.1x at 8 threads
Done
August 4 - August 11
  1. Rework the r.proj band sizing after review feedback that the search approach was hard to follow, replacing it with a footprint grid computed once when the user runs the module
  2. Fix a problem related to low memory slowdown on strongly curved projections and fix it by choosing band height and tile width together from the grid
  3. Verify that the output is correct across all methods, thread counts and memory settings, run the full benchmark matrix, open draft PR #7807
  4. Rebuild the r.geomorphon pytest suite on a DEM that produces all ten landform classes and address review (PR #7785)
  5. Address review on PR #7783
Done
August 12 - August 18
  1. Address remaining review on PR #7785, merged
  2. Rebase PR #7783 over main and add the parallel identity tests
  3. Address any review or questions on #7807 and #7766
  4. Update the wiki page with the final report
Done
Final Week August 19 - August 26
  1. Finish the final report and share it with mentors
  2. Submit the final work product and final evaluation
  3. Respond to review on the open PRs
Done

Reports

Community Bonding Period

During this period I:

  1. Started off with an introductory call with my mentors.
  2. I introduced myself to the community on the OSGeo discourse.
  3. I started this wiki page, created the abstract, project scope, and the timeline.
  4. Did a thread safety audit of the gproj library calls
  5. Studied how the code in readcell.c works
  6. Set up the benchmarking infrastructure and figured out how it works
  7. Looked at the remaining modules and understood how they work
  8. Finalized my dev environment
  9. Confirmed my benchmarks are reproducible
  10. Decided to parallelize r.param.scale first, r.proj second, and then finally r.geomorphon
  11. Worked on identifying methods to parallelize r.param.scale

Week 1

Posted what I did in week 1 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-1-report/153925

Week 2

Posted what I did in week 2 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-2-report/154036

Week 3

Posted what I did in week 3 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-3-report/154478

Week 4 & Week 5

Posted what I did in weeks 4 and 5 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-4-5-report/154479

Week 6

Posted what I did in week 6 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-6-report/154480

Week 7

Posted what I did in week 7 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-7-report/154481

Week 8

Posted what I did in week 8 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-8-report/154560

Week 9

Posted what I did in week 9 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-9-report/154659

Week 10

Posted what I did in week 10 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-10-report/154761

Week 11

Posted what I did in week 11 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-11-report/154833

Week 12

Posted what I did in week 12 in the OSGeo discourse: https://discourse.osgeo.org/t/gsoc-coding-period-week-12-report/154906


Final Report

Abstract

Many GRASS raster modules still do their work on a single core. On a laptop with eight cores that means seven of them sit idle, while only one core does all the work. On a large map, only using one thread means the wait time to receive the output will be long. This project parallelizes three commonly used modules so they can use multiple cores at once, cutting wait time significantly. The modules are r.param.scale, r.geomorphon and r.proj. Each of the three modules is parallelized using OpenMP and a memory bounded approach that gives users the freedom to choose how much memory they want to use. As a result of this project, users can now benefit from a speedup depending on the number of threads they allow and the amount of memory they allocate. The outputs are exactly the same as when the modules ran single threaded, bit for bit.

My Contributions

r.param.scale

r.param.scale is a module that takes an elevation map and computes terrain parameters like slope and curvature at every cell. For each cell in the input map, it looks at the cell’s neighbors and takes their height to find a curved shape that passes through them. The module finds the terrain parameters from this shape.

Each cell's answer depends only on its neighbors, so the map can be split across threads and each thread can work through its share independently. This made the parallelization process simple, but the problem was optimizing memory. The one threaded version of the module loaded the whole map into RAM at once. That was fine for a small map and a single thread, but it doesn’t scale. I rewrote that logic to work through the map in horizontal bands, where each band is sized so it fits inside whatever memory the user allows. The threads then split the rows inside each band. This gives around 4.6x speedup compared to serial at eight threads.

While testing the parallelization, I came across a bug that existed since 2012 within the core GRASS math library. There was a race condition between two shared variables that caused wrong answers to occur when multiple threads were called at the same time. I diagnosed the issue, found a solution, tested it, and added the bug fix within the same PR that parallelized r.param.scale.

Number of threads Speedup
1 1.0x
2 1.91x
4 3.37x
8 4.55x

r.geomorphon

r.geomorphon takes an elevation map as input and labels every cell as a landform type like ridge, valley, slope or flat. It does this by looking outward from each cell in eight different directions and checking if the ground rises or falls along each of the directions. The pattern of ups and downs helps decide what type of landform the cell is.

Almost all of the run time was in that one step, and each cell's answer only depends on the terrain around it. So, I split the output rows across threads. To keep the memory usage limited to how much the user allows, the module now works through the map in horizontal bands. A horizontal band is a certain number of rows clumped together all handled by threads. The size of each band depends on how much memory the user allows. A cell near the top or bottom of a band still needs to see the cells just outside the band. This is because its landform depends on the terrain around it. So each band is read with a few extra rows above and below it. That way every cell sees exactly the same neighbors it would if the whole map were loaded. Not only does this allow the output to be correct, but it also gives around 5.2x speedup at eight threads.

Number of threads Speedup
1 1.0x
2 1.9x
4 3.5x
8 5.3x

r.proj

r.proj is a module that takes a map made in one coordinate system and converts it into another. This is important because many GIS projects need maps with the same projection so that they can take accurate measurements, process data more easily, or align raster data. To build the new map, r.proj takes each cell of the new map, finds the spot in the old map that it corresponds to, and copies the value from there.

Each output cell can be computed on its own, so the rows can be split across threads. So similar to the other 2 modules, the parallelization is simple, but memory is the problem again. Because the two maps are in different projections, one output row does not match one input row. It touches a curved band of input rows, and to compute a set of output rows you have to hold all the input rows they touch. That number changes with the projection and with where you are in the map. On top of that, you have to know where you land before you read the input.

Instead of working this out band by band, the module figures it out once at the start. The module projects the edges of the output map into the input map. Then for each output row, and for each block of columns within that row, it records the input rows that block needs and stores it within a table. That table is used to look things up. After that, sizing a band properly is just looking up the values from the table. It takes as many output rows as it can fit inside the memory the user allows, and if even one full row is too wide it cuts the row into column tiles. If nothing fits at all it uses the original serial code instead. Then each band's input rows are read once and the threads split up the output rows. This makes the speedup range from 2.5x to 5.7x depending on the map and the memory allowed.

One of the many r.proj scenarios: EPSG:4326 to EPSG:3857, memory=300, method=nearest

Number of threads Speedup
1 1.27x
2 2.17x
4 3.54x
8 4.19x

Pytest and Benchmarks

I replaced all three modules' testsuite with a new pytest and added a benchmarking script within the module.

Future work

There is still plenty left to do here:

  • Many raster modules still run single threaded, so they need to be parallelized. The approach used in these three modules transfers directly to any module that works row by row.
  • The progress reporting function in the GRASS library is not thread safe, so parallel modules currently have to work around it. A library level fix would clean that up for every module at once.
  • Several modules still use the old testsuite framework and need to move towards the pytest format.

Conclusion

Three of the more popular GRASS raster modules now run on multiple cores instead of one, and they do it without changing a single value in the output. Users pick how many threads they want to use and how much memory to give them, and the module does the rest.

Along the way I encountered some bugs and was able to fix them. Two of the more prominent bugs were data races in GRASS libraries. One of them sitting there since 2012, and a crash in r.geomorphon on small regions. All three modules now have pytest suites and published benchmarks, so anyone can check the numbers and build on the work.

I want to thank my mentors Anna Petrasova and Huidae Cho for their guidance and reviews throughout the summer, and the GRASS community for being welcoming the whole way.

Log of Pull Requests

Pull Request / Issue Description Status
PR #7440 r.param.scale parallelization and G_ludcmp race fix Merged
PR #7764 Fix a data race on two globals in the projection library Merged
PR #7773 Fix a r.geomorphon crash on regions smaller than the search window Merged
PR #7785 Replace the r.geomorphon testsuite with pytest tests Merged
PR #7627 r.proj parallelization with memory bounded bands Closed in favor of #7807
PR #7807 Simpler r.proj band sizing from a precomputed footprint grid In review
PR #7783 r.geomorphon parallelization Merged
PR #7766 Replace the r.proj method testsuite with pytest tests Merged
Issue #7539 G_ludcmp data race report Addressed within #7440
PR #7185 r.proj proof of concept Closed in favor of #7627
PR #7236 r.param.scale proof of concept Closed in favor of #7440