topic Re: Thinning points based on a time interval in ArcGIS Pro Questions

Thinning points based on a time interval

JeremyDavies — Fri, 10 Apr 2020 17:38:42 GMT

Not sure of the etiquette here, but I am reposting this question on the ArcGIS Pro feed, as I didn't get any responses under the ArcMap heading.

************

Hi-- I want to reduce the number of points in vessel tracklines based on a time interval from a datetimestamp value. I can reduce the number of points based on distance using the 'Delete Identical' tool, but I'd like to be able to reduce the number of points by selecting points based on a user defined interval, like a two or five minute interval, and culling the extraneous points. I have not found any obvious way to do this with the built in tools in ArcMap or Pro, but I'm hoping someone else has worked through this and can suggest a workflow or solution. I am still more comfortable working in ArcMap, but can work in Pro if the solution lies there.

Here is simplified version of what the tabular data look like. The highlighted rows show what records would be selected if I chose four minute interval. It makes sense that some rounding of the times may be needed to standardize the time breaks.

Here is what the desired output might look like.

Thanks!

Re: Thinning points based on a time interval

DanPatterson_Retired — Sat, 11 Dec 2021 08:47:08 GMT

You could do a cumulative sum of the field, take the integer of that and the modulos

A poor example follows... 0 represents the points to keep. Query for those.

a = np.random.randint(0, 10, 20)
s = np.cumsum(a)
m = s % 4

a
array([7, 5, 3, 3, 5, 1, 2, 7, 7, 6, 9, 0, 6, 7, 6, 6, 4, 6, 7, 9])

s
array([  7,  12,  15,  18,  23,  24,  26,  33,  40,  46,  55,  55,  61,
        68,  74,  80,  84,  90,  97, 106], dtype=int32)

m
array([3, 0, 3, 2, 3, 0, 2, 1, 0, 2, 3, 3, 1, 0, 2, 0, 0, 2, 1, 2],
      dtype=int32)

# ---- sample where values are '0'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Thinning points based on a time interval

DanPatterson_Retired — Sat, 11 Dec 2021 08:47:11 GMT

cumulative sum options... python parser.

"""-----------------------------------------
Return the cumulative sum of a field
cumulative(field_name)  # enter into the expression box
"""
old = 0   # include this line
def cumulative(new):
    '''accumulate values'''
    global old
    if old >= 0:
        old = old + new
    else:
        old = new
    return old

#  or

total = 0
def cumsum(in_field):
    global total
    total += in_field
    return total‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Thinning points based on a time interval

JeremyDavies — Fri, 24 Apr 2020 18:52:45 GMT

Thanks for the input, Dan! My computer is in for much needed maintenance, but I'll ponder and try this approach when I am back at it next week.

Re: Thinning points based on a time interval

JoshuaBixby — Sat, 11 Dec 2021 08:47:14 GMT

The following in the interactive Python window will create a new feature layer or table view and make a selection against it:

from arcpy import da, Describe
from arcpy.management import MakeTableView, MakeFeatureLayer
from itertools import groupby

interval = 3 # interval between 1 and 60 minutes
tbl = # path to feature class or table
case_fields = ["Vessel"]
time_field = "DTS_DateField"

oid_field = Describe(tbl).OIDFieldName
fields = case_fields + [time_field, oid_field]
sql_orderby = "ORDER BY {}, {}".format(", ".join(case_fields), time_field)
time_group = lambda dt: dt.replace(minute=dt.minute//interval*interval, second=0, microsecond=0)

oids = []
with da.SearchCursor(tbl, fields, sql_clause=(None, sql_orderby)) as cur:
    key_func = lambda row: row[:-2] + (time_group(row[-2]),)
    for key, group in groupby(cur, key_func):
        oids.append(next(group)[-1])
        

if Describe(tbl).datasetType == 'FeatureClass':
    lyr = MakeFeatureLayer(tbl, "lyr").getOutput(0)
elif Describe(tbl).datasetType == 'Table':
    lyr = MakeTableView(tbl, "lyr").getOutput(0)
else:
    exit()

lyr.setSelectionSet(oidList=oids, method="NEW")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The above selects the first record in a 3-minute interval. You can pick whatever minute interval you want. It is also possibly to get fancier with record selection within the time interval, but selecting the first is the simplest.

Re: Thinning points based on a time interval

JeremyDavies — Fri, 01 May 2020 19:26:57 GMT

Thanks Joshua! I just tested it out and it works great for thinning points at whatever minute interval I choose. This could be a game changer for me, assuming I can handle some of the massive datasets I have.

Re: Thinning points based on a time interval

JeremyDavies — Sat, 11 Dec 2021 08:47:17 GMT

Hey Joshua-- I'm not sure if you'll be willing to chime in on this, but strangely I am not able to run this script today using the same data, project and syntax as last week. I am not well versed in python (I plan to take esri's online course soon!), so I'm guessing the solution is simple and just reflects my ignorance. If you want, I can send you a small sample dataset and you can see what I am working with. But I understand if you don't have the time or desire to offer more me free help!

Below is the python script that I successfully last week, where the input point layer (tbl = “test_points”) lives in the project's default database (C:\GIS\geodata\USA\Vessel Presence Study\AIS Data\MXAK Data\MXAK Sample Data.gdb).

from arcpy import da, Describe 
from arcpy.management import MakeTableView, MakeFeatureLayer 
from itertools import groupby 
  
interval = 10 # interval between 1 and 60 minutes 
tbl = “test_points” 
case_fields = ["MMSI"] 
time_field = "dtg_test" 
  
oid_field = Describe(tbl).OIDFieldName 
fields = case_fields + [time_field, oid_field] 
sql_orderby = "ORDER BY {}, {}".format(", ".join(case_fields), time_field) 
time_group = lambda dt: dt.replace(minute=dt.minute//interval*interval, second=0, microsecond=0) 
  
oids = [] 
with da.SearchCursor(tbl, fields, sql_clause=(None, sql_orderby)) as cur: 
    key_func = lambda row: row[:-2] + (time_group(row[-2]),) 
    for key, group in groupby(cur, key_func): 
        oids.append(next(group)[-1]) 
         
 
if Describe(tbl).datasetType == 'FeatureClass': 
    lyr = MakeFeatureLayer(tbl, "lyr").getOutput(0) 
elif Describe(tbl).datasetType == 'Table': 
    lyr = MakeTableView(tbl, "lyr").getOutput(0) 
else: 
    exit() 
  
lyr.setSelectionSet(oidList=oids, method="NEW") ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I get the following syntax error, even though this ran before:

One related question, assuming I can get the script to runs: can I write out the full path for 'tbl', as in

'tbl = C:\GIS\geodata\USA\Vessel Presence Study\AIS Data\MXAK Data\MXAK Sample Data.gdb\test_points' , where 'test_points' is the point data layer I want to thin? I haven't had any luck with that syntax.

Thanks for any input you may have! And if you want to me to send a small dataset so you can see what I'm working with, you can email me at jdavies@oceanconservancy.org, unless there is a way for me to share data with via geonet. And if you want me to work this out on my own, just say the word!

Re: Thinning points based on a time interval

JoshuaBixby — Thu, 07 May 2020 15:18:33 GMT

First, let's get the SyntaxError error resolved, then you can follow-up with another question. You are getting the SyntaxError because of your quotes. If you look closely at your quotes, you are using left double quote and right double quote, which are not valid. You need to use the regular ASCII double quotes.

Re: Thinning points based on a time interval

JeremyDavies — Thu, 07 May 2020 20:23:33 GMT

Yep- that did the trick! I switched to using Notepad++ for editing, and now it runs fine. I suspect this was the quotes were the problem with using the full path as well. Thanks!

Re: Thinning points based on a time interval

JeremyDavies — Fri, 20 Nov 2020 18:54:32 GMT

Hi Joshua (or any other python gurus who see this)-- this python script has worked great for me for months, but for some reason I keep getting the following error when I run it today, even using the exact same data and syntax that ran perfectly yesterday:

Traceback (most recent call last):
File "<string>", line 29, in <module>
AttributeError: 'str' object has no attribute 'setSelectionSet'

I've tried restarting my project as well as creating a new project, but to no avail.

As a reminder, here is the code again, and I've attached a simplified gdb that includes a stripped down version of some of the point data. If have you have any insights, I'd really appreciate it!

Thanks-- JD

from arcpy import da, Describe from arcpy.management import MakeTableView, MakeFeatureLayer from itertools import groupby interval = 10 # interval between 1 and 60 minutes tbl = "MXAK_Arctic_2019_11_AllPoints" case_fields = ["MMSI"] time_field = "Base_station_time_stamp" oid_field = Describe(tbl).OIDFieldName fields = case_fields + [time_field, oid_field] sql_orderby = "ORDER BY {}, {}".format(", ".join(case_fields), time_field) time_group = lambda dt: dt.replace(minute=dt.minute//interval*interval, second=0, microsecond=0) oids = [] with da.SearchCursor(tbl, fields, sql_clause=(None, sql_orderby)) as cur: key_func = lambda row: row[:-2] + (time_group(row[-2]),) for key, group in groupby(cur, key_func): oids.append(next(group)[-1]) if Describe(tbl).datasetType == 'FeatureClass': lyr = MakeFeatureLayer(tbl, "lyr").getOutput(0) elif Describe(tbl).datasetType == 'Table': lyr = MakeTableView(tbl, "lyr").getOutput(0) else: exit() lyr.setSelectionSet(oidList=oids, method="NEW")

Re: Thinning points based on a time interval

TanuHoque — Wed, 31 Mar 2021 03:35:04 GMT

I somehow stumbled on this post from late last year. I thought I should provide another solution to this problem for folks from the future 🙂

Pros: the benefit of this approach is that this dynamic. When new data comes in, it will continue to work without you doing anything.
Cons: you must use enterprise database. FileGDB won't work.
finally the sql that I provided below is not a very elegant one. I put that together quickly and also helps you understand what is going on in this nested fashion.

Let's take a look at the output first:

Black dots represents the original vessel's positions in the screenshot below. Red circles, off a different layer, show positions with some temporal gaps. Imagine you group the original locations (black dots) into bins/buckets of 5 min interval. Then pick up first ones (red circles) from each group. That is what I did here.

How to:

first you must copy data in an enterprise database (in my case I copied them in a SQL Server)
i used the data attached in the previous comment
click Query Layer from Add Data dropdown in Map tab in ArcGIS Pro
choose your database connection
provide a name
paste the following SQL (this is for SQL Server, and should work with other database flavors too)

select objectid, datetimestamp, mmsi, latitude, longitude, shape from ( select *, rank() over (partition by elapsed_time_group order by datetimestamp) as rnk from ( select *, DATEDIFF(second, FIRST_VALUE(datetimestamp) over(order by datetimestamp), datetimestamp) / 300 as elapsed_time_group from VESSELTRACK_MMSIPOINTS ) a ) b where b.rnk = 1

click Validate
click Next
make sure the correct field is chosen as the Unique Identifier Fields(s) aka the object-id field
make sure the Geometry type and spatial reference is picked correctly. Otherwise choose one.
click Finish
a query layer gets added to your map

Problem with this SQL is that it produces a result where a few positions are still very close together (see the screenshots below). That is because positions are not recorded at a regular frequency.

To weed out those points too, I have added an additional filter in the SQL - see below.

select objectid, datetimestamp, mmsi, latitude, longitude, shape from ( select *, datediff(minute, lag(datetimestamp) over (order by datetimestamp), datetimestamp) as time_diff from ( select objectid, datetimestamp, mmsi, latitude, longitude, shape from ( select *, rank() over (partition by elapsed_time_group order by datetimestamp) as rnk from ( select *, DATEDIFF(second, FIRST_VALUE(datetimestamp) over(order by datetimestamp), datetimestamp) / 300 as elapsed_time_group from VESSELTRACK_MMSIPOINTS ) a ) b where b.rnk = 1 ) c ) d where d.time_diff > 2 or d.time_diff is null

here is the output where green circles represents a cleaner version of the previous one.

Hope some folks might find this helpful in the future 🙂