Computing 3D seismic fold is quite straight forward, for those unfamiliar with the concept, below is a brief explanation.

Background information (Seismic 101) :

Seismic Fold:

Basically, I have (several hundred thousands of) Receiver Points (RP) and Source Points (SP). Each Source, when activated, will send a signal that will be picked up by (several thousands of) Receivers.
For every pair (SP, RP) I will compute the Midpoint, i.e. the coordinates of the middle of the segment (the ray) joining the 2 points and find out in which bin this Midpoint falls.
The bins grid is a regular grid (usually orthogonal/rectangular) with cell dimensions being RP_i/2 and SP_i/2 where _i denotes the interval between points (these intervals are fixed as they are taken on the theoretical RP/SP grid).
The Fold Map will then simply be a heatmap with the values being the number of Midpoints falling into each bin of the bins grid.
Ideally, I'll need to be able to filter the Midpoints based on the Offset length. The offset being the length of the ray (i.e. the distance between SP and RP for the pair considered).
A typical fold map may look like this.

SPS (Shell Processing Support) files:

The information source is usually in the form of 3 SPS files (see format definition here). I am able to load the 2 Point files (R and S) into QGIS so I get all Receivers and Sources positions.
The Relation (X) file gives the list of all RP listening to each SP.

The problem:

The problem lies mainly in the numbers. I have to compute several million rays so every nano-second will count here!
The X file lists the points by name (or rather, number) so I'll have to lookup the coordinates of these points every single time as I'm reading through the X file.
Is there an efficient way to load my R and S files in memory (ideally caching only a portion of these files at a time) and how should I hold the information of the computed Midpoints ? In memory ? In a numpy array ? Is there a QGIS object that would be particularly suitable for this kind of extensive computation ?

  • I have been working on this for the past year or so sporadically. I had originally written code within a Visual Foxpro Database, and decided to translate the code into either C++, or Python. I get sidetracked by work quite frequently so I only have the mid-point (CDP) coordinate calc routines finished. This was the simple part. I believe from my experience it will be faster to create routines that operate outside of QGIS, creating WKT files of the grid, CDP locations, and relations files, then creating watch files within QGIS for the updates to appear on screen.
    – jbgramm
    Commented Dec 2, 2018 at 20:26
  • Other things to consider are Hexagonal bins, rectangular rather than square bins, overlapping bin boundaries (Decrease the bin size by 0.0001 units of measurement to eliminate a CDP falling in two bins) Simplify the calculations so all of the trig is taken care of early, then use simple math based upon the trig constants for the rest of the calculations. There are Hexagonal, and regular grid generation utilities available within the processing toolbox which may give you a foundation to work with.
    – jbgramm
    Commented Dec 2, 2018 at 20:36
  • I'm hoping to do as little computation as possible, leaving numpy doing the binning for me if possible. The main bottleneck I see from where I am today (not very far) is calling the coordinates of each point.
    – YeO
    Commented Dec 2, 2018 at 21:45
  • Have a look here: agilescientific.com/blog/2015/1/8/it-goes-in-the-bin . Do not miss the first comment as it's using a more optimised geopandas tool. Commented May 30, 2019 at 8:49
  • While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review Commented May 30, 2019 at 9:14

1 Answer 1


This is too long for a comment, but it is not the answer to your issues.

Unfortunately, there is no avoiding that bottleneck.

You could theoretically set limits to the calculation ranges by point numbers (Line & Station, or Track & Bin) if you know the spacing, and make general assumptions with regards to the distance from a receiver point to a source point based upon station interval, and line spacing

You obviously know the following, but I will expound on it anyway.

You are going to have two loops, or do while statements.

One will fetch each receiver point (R File), storing the coordinate values, and other pertinent information.

The second loop will fetch each source point (S File), calculate the inverse or ray distance, and also take the x, and y coordinates of both points, sum the x values and the y values then divide each sum by 2 to derive the mid point coordinates.

I would also take the inverse or ray distance, and divide this by two, this will give you the depth of the record, and I would calculate the azimuth of the ray.

Then write all of these values to another file. (The relations or X file).

So in the third file I would write the Receiver Point Number, Source Point Number, MidpointX, MidpointY, Radial_distance, Half_Distance, Azimuth.

You could add other info if you needed to such as Receiver Coordinates, and Source Coordinates, Preplot or Postplot Designators and so on.

Then after each source point has been looped through, you return to the receiver point file, fetch the next one (first loop), then repeat the process with the source points (second loop) until all of the points in both files are looped through.

It is very calculation intensive, but the math is simple for the most part.

Optimization of the loops is one of the most important things here. On a seismic grid with 8000 receiver points, and 8000 source points, you will have 64,000,000 mid point locations, and in all honesty, that is not a very intense grid compared to some I have worked on.

The rest of this is not meant to discourage you, but I see the road you are taking, and it will progress to the following.

There is a lot more fun to be had when you start to analyze azimuthal distribution, and create spider graphs for each bin. Then querying specific depth ranges for fold.

Other fun will be the blending of surveyed locations with unsurveyed locations for design purposes and to avoid the striping of the data.

It can then turn into something even more mind boggling when you start turning off specific points within the grid, and have to regenerate the counts for the mid point locations within specific bins based upon the loss of these points.

Have as much fun with this as you can. Take small steps, and work on optimizing the loops first.

So in a follow up to this, out of morbid curiosity, insanity, and stubbornness, I decided to try this using only the core tools, and processing tools contained within QGIS 3.4.3.

I will not go into full step by step details, but I started by using the regular points generator, and built a grid of points on a 220 ft. by 220 ft. grid for both receiver points, and source points. I then rotated both grids using the processing toolbox rotate function.

After that I built a polygon grid of 110 ft by 110 ft, and then translated and rotated the grid so the points were centered within the polygons.

I then queried out the desired points from the receiver and source point files.

Picture of the Polygon Grid

Grid and Preplots

Then I ran the Distance Matrix utility from the Vector Menu (it is also part of the Processing Toolbox)

This created a record for each point on the Receiver Point Grid related to each point on the Source Point Grid. In all there are 3,128,160 combinations. The distance matrix file gives you an input and target point id, and the distance between the two points, which in many cases was found to be erroneous unless I am misunderstanding what is happening there.

I then joined the receiver point file with the distance matrix file exporting the coordinates of the receiver points associated with the file.

I repeated the above process using, this time joining the source point file with exported file.

From there I added a mid point x and y field, as well as a mid point distance field, and updated those fields using the SQL functions. This was a very time consuming process with frequent crashes along the way. I attempted this using both shape files, and spatialite files. The spatialite files seemed to be more stable at least initially. (After each crash, I rebooted the computer.)

I took the file with the calculated midpoint data, and exported a CSV file of containing the Receiver Point ID, Source Point ID, Mid Point X, Mid Point Y, and the Mid Point Distance.

I then opened Mid Point CSV file, and exported it as a Spatialite file. The point count was correct, however there seems to be an error in the coordinates as some of the points fall outside of the boundary of the Bin Grid.

Mid Points on Bin Grid

I have not tracked this error down as of yet by comparing the joined receiver and source grid coordinates with the original coordinates, that is next on my list.

Out of curiosity I decided to calculate a full fold map with what I came up with. I knew it would be wrong, but I wanted to see how long it would take counting the points in polygons. (Using the first run shape files)

It took a little over eight hours using the Points in Polygons routine.

I am not going to run it again using the spatialite files until I solve the rotational issues which is apparent in both the shape files of the first run, and the spatialite files of the second run.

I am still having problems wrapping my head around using Python and looping through both the receiver and source point files staying with a single point in the receiver file while calculating mid points and distances to each source point from that receiver point, then moving to the next receiver point and repeating this process until all points are calculated.

The rest of the math is easy, but the looping thing has given me fits. All I want to write out to a CSV file is the Receiver Point ID, Source Point ID, Mid Point X & Y coordinates, and the Distance or Depth. That is the easy part for me so far.

Given what I have experienced with the above experiment, I hope you have a truly monstrous computer, huge hard drives, and a ridiculous amount of RAM, along with an inordinate amount of patience.

  • I'm thinking down that road too... but for now, it would be nice to be able to compute even the basic full fold (all offsets). From a theoretical grid, we can probably save a lot of computation by using grid locations, but I'd like to be able to compute from a postplot also and use actual locations (you cannot assume grid locations then). Also, the project I'm looking at right now has about 800,000 SP and 6,656 RP for each so we're talking 5.3 billion midpoints... I'll keep studying what numpy has to offer as I'm sure there are a few functions that can handle most of the heavy lifting.
    – YeO
    Commented Dec 3, 2018 at 10:00
  • With that density of source points, are these in a source point array where you can calculate a center of the group rather than using each individual source point location? If you are able, that is the way I would approach this depending on the nominal fold desired.
    – jbgramm
    Commented Dec 3, 2018 at 15:09
  • Nope, these are single sources (groups of 1) so there is no shortcut here...
    – YeO
    Commented Dec 3, 2018 at 15:16
  • Ouch, or wow or something. There is probably a simple algorithm which would give you preplot full fold based upon grid number separation between the receiver and source points but that will only give theoretical, or design fold, and not real world fold. My lack of familiarity with Numpy puts me at a disadvantage here. The brute force routines that I wrote to calculate midpoints took about 5 minutes for a 10000 by 10000 grid, but it wrote a lot of data beyond the midpoint locations, and it was on a slow computer. I would not want to calculate the fold on your grid with that computer.
    – jbgramm
    Commented Dec 3, 2018 at 16:56
  • I believe preplot full fold should be easy as it is for the main part applying the same pattern over and over and just identifying edge cases where not all offsets are present. That should be manageable. But real world data is where it gets interesting...
    – YeO
    Commented Dec 3, 2018 at 23:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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