Is there a way I could program my own script without too much effort?
The most effortless way is running the tools manually, then going to the Geoprocessing history and copying the commands into a Python script:
I found another way that might help you:
- Calculate the distance of your measured points to the transmitter.
- Generate buffers around the points, based on that distance.
- Get the intersection points of all the buffers
- In a perfect world without measuring uncertainties and without signal absorption in walls, all circles around your points with a radius of the calculated distance would intersect in a single point, the transmitter (and in many other points, but they would all have that common point).
- In our more imperfect world, the highest density of intersections should still be close to the transmitter.
- Create a heat map of the intersection points.
- Analyze that heat map like in my original answer.
Just to show how you could create your own script, I'll run the tools manually and then copy/paste the Python commands from the history here. Note that I'm using memory paths for intermediate results. This means they are only stored in RAM and thus will be gone when I close ArcGIS. The advantage is that this way is much faster than saving on disk.
My input data, area is ~ 250m x 250m, signal strength is calculated using the formulas from the links below plus a little noise.
Calculate the distance of your measured points to the transmitter. In this thread and this wiki page, you can find formulas to convert signal strength to distance to the transmitter.
arcpy.management.CalculateField("points", "Distance", "dist(!DoubleField!)", "PYTHON3", """def dist(fspl):
k = -27.55
f = 2413
pot = (fspl - k - 20 * math.log10(f)) / 20
return math.pow(10, pot)""", "TEXT", "NO_ENFORCE_DOMAINS")
Create buffers around the points based on the calculated distance.
arcpy.analysis.Buffer(r"memory\points", r"memory\circles", "Distance", "FULL", "ROUND", "NONE", None, "PLANAR")
Intersect the buffers with themselves, specify point output.
arcpy.analysis.Intersect(r"memory\circles #", r"memory\circle_intersections", "ALL", None, "POINT")
Create a heat map of these intersections (Point Density). Play with the Output cell size and Neighbourhood parameters, they massively influence the result.
out_raster = arcpy.sa.PointDensity("circle_intersections", "NONE", 1, "Circle 5 CELL", "SQUARE_METERS"); out_raster.save(r"memory\heatmap")
The area with the most intersections of the buffers should be near the transmitter. Create Contour lines to find this area.
arcpy.ddd.Contour("heatmap", r"memory\contours", 1, 0, 1, "CONTOUR", None)
That's pretty busy. Let's zoom in.
Slightly over 5m off, not too shabby!
Let's create an ouput feature class (saved on disc), showing the most probable transmitter locations.
# get the raster values as points
arcpy.conversion.RasterToPoint("heatmap", r"memory\raster_values", "Value")
# delete all points where grid_value != max(grid_value)
# it would be easier to use the arcpy.da.SearchCursor and UpdateCursor for this (or do it manually in the attribute table)
# but we're pretending we don't know Python and want to create a script, so let's do it with tools.
arcpy.management.CalculateField("raster_values", "IsMax", """if($feature.grid_code == 0) {
return 0
}
var max_grid_code = Max($featureset, "grid_code")
return $feature.grid_code == max_grid_code""", "ARCADE", '', "TEXT", "NO_ENFORCE_DOMAINS")
# 400k points take a really long time...
# 1/2 hour later: well, I'm going to use the Cursors...
# this is how you would procee with tools:
arcpy.management.SelectLayerByAttribute("raster_values", "NEW_SELECTION", "IsMax = '1'")
arcpy.management.DeleteFeatures("raster_values")
# And this is using cursors
grid_codes = [row[0] for row in arcpy.da.SearchCursor("memory/raster_values", ["grid_code"], "grid_code > 0")]
max_grid_code = max(grid_codes)
with arcpy.da.UpdateCursor("memory/raster_values", ["grid_code"], f"grid_code < {max_grid_code}") as cursor:
for row in cursor:
cursor.deleteRow()
# that took like 10 seconds...
# save the result on disc
save_path = "TransmitterLocations" # Save in the project's default gdb
arcpy.management.CopyFeatures(r"memory\raster_values", save_path)
Have a great day!
Johannes