# Number of Datapoints / Area

**URL:** <https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765>\
**Category:** Q&A\
**Created:** [August 12, 2020, 8:18am UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765 "2020-08-12T08:18:20Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 12, 2020, 8:18am UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/1 "2020-08-12T08:18:20Z")

</div>

Hi All,

I have generated a global dataset of mountain peaks (arbitrarily located (x,y,z) triples) in a geographic coordinate systems (WGS84 - long/ lat). Now I try generating a global grid by computing the number of peaks / area. It seems that this should be a simple task for GMT but I cannot find a working approach (sorry for the trivial question). I would require counting the number of peaks within a moving circular window with a given radius (in km).

Any ideas how to perform this task?

Cheers Jörg

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 12, 2020, 6:44pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/2 "2020-08-12T18:44:16Z")

</div>

There is no specific module that does exactly what you want. You can do it with a square box in lon/lat via blockmean or xyz2grd. Or, you can write a script that uses gmtselect to find the number of your xyz triplets taht are inside the circle centered on each desired grid node in a loop. Should not be too hard: Pick your desired global grid parameters (increment in lon,lat), use grdmath to make an empty grid, dump it to ascii, then loop over the records and make a single-record point file that you use in gmt select -C to find the points inside that radius and then add lon lat count to a growing output table that you then xyz2grd make into a grid in the end.

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 12, 2020, 7:12pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/3 "2020-08-12T19:12:51Z")

</div>

Feel free to comment on my feature request at [https://github.com/GenericMappingTools/gmt/issues/3942](https://github.com/GenericMappingTools/gmt/issues/3942) if I missed a functionality.

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 12, 2020, 8:03pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/4 "2020-08-12T20:03:14Z")

</div>

Dear Paul,

Thank you for showing a working solution for this problem. I will give it a try and write a little script following your recipe.

The feature request perfectly describes the missing functionality. I was thinking that nearneighbor could do the job but there is no option for writing the sum the of data points within the search radius and as far as I understand it just takes the nearest neighbor of each sector and not all data points within the search radius. Anyhow, gmtcount sounds great!

Cheers  
Jörg

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 13, 2020, 3:33pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/5 "2020-08-13T15:33:44Z")

</div>

Hi Paul, colleagues,

I have translated Paul’s recipe into a small bash script. As expected the bash script has a bad performance and is not suitable for a large scale / high-resolution analysis. Anyhow, maybe this workaround is useful for somebody else. I think slicing large domains into smaller tiles by a concurrent reduction of data points should increase the speed and may allow a simple form of parallelization. Any further ideas to improve the performance?

Here is the bash script

#Recipe of Paul  
#-----------------------------------------------------------------------------------------------------  
#Decide on your grid specifications and make a dummy empty grid with grdmath  
#Dump the grid with grd2xyz to ascii  
#Write a script that loops over the records in this list of nodes  
#For each record, dump it to a point.txt file  
#Run gmtselect -Cpoint.txt +dradius on the whole dataset, then count the number of records found  
#Add the node location and the count to a growing output table  
#Run xyz2grd on that table to make the final grid

function NumberOfDataPoints()  
{

grdmath $reg $res -r 0 = dummy.grd  
grd2xyz dummy.grd | awk ‘{print $1,$2}’ \> Nodes

echo lon lat NoP \> Density.txt  
while read line; do  
echo “$line” \> point.txt  
NumberOfPeaks=$(gmtselect peaks -Cpoint.txt+d"$radius"km -fg | wc -l)  
echo “$line” $NumberOfPeaks \>\> Density.txt  
echo “$line” $NumberOfPeaks  
done \< Nodes

xyz2grd Density.txt $reg $res -r -h1 -GDensity.grd

#clean up  
rm dummy.grd Nodes point.txt  
}

#Define Region of interest  
W=0  
E=20  
S=40  
N=50

reg=-R"$W"/"$E"/"$S"/"$N"

#Resolution of the resulting grid  
res=-I30m

#Name of the input dataset  
data=peaks

#Radius of the moving Window in [km]  
radius=100

#Generate Subset of Input Data: Region of Interest + Buffer larger the Radius of gmtselect  
#Buffer in degrees  
Buffer=1

Win=$(echo “$W - $Buffer” | bc -l)  
Ein=$(echo “$E + $Buffer” | bc -l)  
Sin=$(echo “$S - $Buffer” | bc -l)  
Nin=$(echo “$N + $Buffer” | bc -l)  
regInputData=-R"$Win"/"$Ein"/"$Sin"/"$Nin"

#Preprocess the point dataset --\> Col 1 : long; Col 2 : lat  
cat Peaks.P500.Sorted | awk ‘BEGIN{FS=","}{print $5,$4,$6,$7}’ | gmtselect $regInputData \> $data

#Call function to count data within the defined radius  
NumberOfDataPoints

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 13, 2020, 6:16pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/6 "2020-08-13T18:16:43Z")

</div>

FYI, see [https://github.com/GenericMappingTools/gmt/pull/3950](https://github.com/GenericMappingTools/gmt/pull/3950)

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 13, 2020, 7:24pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/7 "2020-08-13T19:24:15Z")

</div>

Wow – this was fast! I’m keen to try it 🙂

Thank you!

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 14, 2020, 7:38am UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/8 "2020-08-14T07:38:14Z")

</div>

Hi all,

I just saw the discussion about the new module on github and wondering if the code is already part of the developement branch of GMT? I’ve just upgraded to “GMT - The Generic Mapping Tools, Version 6.2.0\_fff59e9\_2020.08.13” but can not find a module “gmtcount” or “gmtstats”.

cheers  
Jörg

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 14, 2020, 7:51am UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/9 "2020-08-14T07:51:56Z")

</div>

It is not merged into master but is presently in a branch called gmtcount, which you could check out and build. Likely to change but not tonight.

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 14, 2020, 8:07am UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/10 "2020-08-14T08:07:49Z")

</div>

Ok - thank you for the fast response. It seems that the branch gmtcount requires a username and a password.  
Anyhow - I think I can wait until the module is part of the master branch.

Thank you for all the effort developing GMT over so many years!

-J

---

<div class="post-metadata">

**Author:** ![jrobl](https://avatars.discourse-cdn.com/v4/letter/j/c68b51/32.png) [@jrobl](https://forum.generic-mapping-tools.org/u/jrobl)\
**Post date:** [August 14, 2020, 12:00pm UTC](https://forum.generic-mapping-tools.org/t/number-of-datapoints-area/765/11 "2020-08-14T12:00:26Z")

</div>

> [@jrobl](#):
>
> It seems that the branch gmtcount requires a username and a password.

Ooops - Eventually I could download the branch (used a wrong adress … ). The module is doing exaclty what I was seeking for and it is really fast 🙂. I love this tool.

Thank you!!!  
Jörg

 ![Density](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/1X/9ad1dbbcafce72ddc2eae6818cc7c6d391971397.png)