Select to view content in your preferred language

How can I improve the accuracy, working in Python?

1049
5
11-10-2022 06:28 AM
Labels (1)
TomGeo
by
Occasional Contributor III

I am working on a little task, where I create some points, based on coordinates. Every point is then used to create a couple of polylines, where each point is supposed to be the start point for a set of polylines. The endpoints for the polylines are determined by angle and distance.

So, the outcome should be a point, in the center of some polylines, looking like sunrays.

I create the points and use the created points thereafter as point of origin for my polylines. Only problem, the created point is not where the polylines 'meet' and the polylines do not meet at the same point either.

I wonder how that can be, if the point of origin does not change!

 

import arcpy
import arcgis
from arcgis.gis import GIS
from arcgis.geometry import Point, Polyline, Geometry

bearing_list = [(345,165),(315,135),(285,105),(255,75),(225,45),(195,15)]
distance_2 = 2000
in_sr = 4326
out_sr = 3006

# create featureclass and populate with points from coordinates
points_fc = arcpy.CreateFeatureclass_management(out_name='Sites', geometry_type='POINT', spatial_reference=out_sr)
arcpy.AddField_management(points_fc, 'SiteName', 'TEXT')
point_insert_cursor = arcpy.da.InsertCursor(points_fc, ('SiteName', 'SHAPE@XY'))
pnt_list = []

for item in point_dict:
    y, x = point_dict[item]
    pnt = Point({
        'x': x,
        'y': y,
        'spatialReference':{
            'wkid': in_sr, 
            'latestWkid': in_sr
        }
    })
    out_pnt = pnt.project_as(out_sr)
    pnt_list.append([item, out_pnt])
    row_point = (item, out_pnt.centroid)
    point_insert_cursor.insertRow(row_point)        
        
del point_insert_cursor

# iterate over points, create temporary line featureclass, 
temp_lines = arcpy.CreateFeatureclass_management(out_name='temp_lines', geometry_type='POLYLINE', spatial_reference=out_sr)
arcpy.AddField_management(temp_lines, 'SiteName', 'TEXT')
line_insert_cursor = arcpy.InsertCursor(temp_lines)
line_length = distance_2 + 5

for pnt in pnt_list:
    for bearing_set in bearing_list:
        pnt_array = arcpy.Array()
        x = pnt[1].first_point.x
        y = pnt[1].first_point.y
        pg = arcpy.PointGeometry(arcpy.Point(x, y), arcpy.SpatialReference(out_sr))
        start = pg.pointFromAngleAndDistance(bearing_set[0], line_length)
        start_pnt = arcpy.Point(start.firstPoint.X, start.firstPoint.Y)
        
        end = pg.pointFromAngleAndDistance(bearing_set[1], line_length)
        end_pnt = arcpy.Point(end.firstPoint.X, end.firstPoint.Y)
        
        pnt_array.append(start_pnt)
        pnt_array.append(end_pnt)
        
        line_feature = line_insert_cursor.newRow()
        line_feature.shape = arcpy.Polyline(pnt_array)
        line_feature.setValue("SiteName", pnt[0])
        
        line_insert_cursor.insertRow(line_feature)

del line_insert_cursor

 

 

The result for one of the points: the point is 10 centimeters away from where the lines 'meet'.

TomGeo_0-1668090345117.png

 

And it looks like this at the 'center' of the polylines

TomGeo_1-1668090479919.png

 

- We are living in the 21st century.
GIS moved on and nobody needs a format consisting out of at least three files! No, nobody needs shapefiles, not even for the sake of an exchange format. Folks, use GeoPackage to exchange data with other GIS!
Tags (3)
0 Kudos
5 Replies
DanPatterson
MVP Esteemed Contributor

from 

Polyline—ArcGIS Pro | Documentation

You should specify a coordinate system when creating the polyline

Polyline (inputs, {spatial_reference}, {has_z}, {has_m})

And you may need to use a geodesic densify on the polyline if you are trying to get all the lines to cross at the same spot

In the Methods section 

densify (type, distance, {deviation})

 


... sort of retired...
0 Kudos
TomGeo
by
Occasional Contributor III

Hi @DanPatterson, thanks for the advice about the coordinate system when creating the geometry. It didn’t change anything on the outcome though. 🙂

Neither did density, and the only thing that solved the issue with the point of origin was to create all lines from the point of origin, leaving me with 12 instead of 6 lines.

I assume the culprit is actually the precision of the points, calculated by pointFromAngleAndDistance. Meaning calculating points on a circle with a distance of 2000 meters, and angles of 345 and 165 degrees will give you points, not marking a precise opposite of each other on a circle.
Once you connect these two points by a line, the line has good chances not to cross your point of origin.

For my purposes it is not important that the point where all polylines meet is spot on on top of the point used for the calculations (even though that’s what I would expect). But it is important that all polylines have a common cross point.
Later in the process these polylines are used to segment a buffer around my calculation point into 12 equal circle segments. As long as the polylines do not have a common crossing point, I will create more than 12 segments, where an unknown number of them will simply be slivers, and in the end the 12 expected circle segment differ in their area size.

- We are living in the 21st century.
GIS moved on and nobody needs a format consisting out of at least three files! No, nobody needs shapefiles, not even for the sake of an exchange format. Folks, use GeoPackage to exchange data with other GIS!
0 Kudos
DanPatterson
MVP Esteemed Contributor

understand... I usually ensure that the coordinates are in a defined projected coordinate system to ensure double precision and I generate sectors, rings and circles from a central point, angle and radius(radii) to get what I want.  I blogged about this a number of years ago.

Circles, sectors, rings, buffers, and the n-gons - Esri Community


... sort of retired...
TomGeo
by
Occasional Contributor III

Certainly bookmarked your blog, and I am pretty sure I will steal some of your examples. 😉

- We are living in the 21st century.
GIS moved on and nobody needs a format consisting out of at least three files! No, nobody needs shapefiles, not even for the sake of an exchange format. Folks, use GeoPackage to exchange data with other GIS!
0 Kudos
RogerDunnGIS
Occasional Contributor II

If the points and lines are in the same geodatabase and in the same geodataset, then they likely have the same XY tolerance, so that's good.  Correct me if I'm wrong, community, but what I understand the XY tolerance does is define how many integers fit into a measurement of 1 map unit.  Databases prefer to store and query ints, so the XY tolerance is what breaks those integers down into floats for what we see on our maps, attribute tables, JSON geometries, etc.  For example, if a feature class is created with an XY tolerance of 0.00328083333333333 in Foot_US, that means that a true database value of 1000 will come back to the software and the user as 3.28083333333333 US Feet.  A double/float value of 612,718.82934 will be stored in the database as an int right around 186,757,073.  Am I right, everyone?

If your points and fireworks aren't in the same geodatabase or geodataset, then maybe the difference in XY tolerance is the problem.

Earth geometry doesn't always work out like planar geometry does.  It still bothers me that terms like x and y are used in GIS when, in reality, the y axes converge to points after several thousand miles.  And it's still a mystery what the earth's actual shape is, what with all the changes that occur to it with seismic and earthquakes, constantly moving tectonic plates, and an ocean that crawls up the shore and recedes at night.  As much as I want my GIS to be like graph paper, it just isn't.

0 Kudos