Tagging links with screenlines ids

[malamillo]

Hello folks,

I'm wondering if there's a tool or python code that would allow me to automate tagging links on a road network by doing a spatial join between the auto links or transit segments and a screenlines GIS layer provided as a shape file. The tag value would be stored as an extra attribute in the links or transit segments.

I know this can be done manually using the recent versions of EMME, however I was wondering if someone has cracked a python code that may automate this task so we don't have to start from scratch.

I will welcome your thoughts. Kind regards,

P.s, Thanks for all the great work developing and maintaining the TMG toolbox.

Mauricio Alamillo Planner Systems Analysis and Forecasting Office Ministry of Transportation Policy and Planning Division Transportation Planning Branch

777 Bay St, 7th Floor, Suite 700 Toronto, ON M7A 2J8 Tel. 416-585-6016 Fax. 416-585-7324

[PeterKucirek]

Hi Mauricio,

It's been a few years since I've used it, but there's a tool called "Load Attribute From Polygon" that ought to do the trick.

P.S. I recommend you edit your post to NOT contain your full contact info as this is a public website.

[malamillo]

Thanks Peter, I'll look into that.

My contact info is already public through the OPS directory. I do appreciate your concern.

[JamesVaughan]

@malamillo Was the "Load Attribute From Polygon" able to solve your use case?

[malamillo]

Hi James.

INRO suggested using the Data Tables API released in 4.3.0. I haven't tried it yet, but it seems to have all the tools needed to carry out spatial join analysis. I'll let you know the results once I get to do some progress on my end.

Regards.

Mauricio

[beykaei]

Hi Mauricio There are some libraries in python such as "PySAL" "fiona" can do spatial analysis including spatial join. Writing your own code is not that much hard using such libraries.

This thread at Stackexchange may help you. https://gis.stackexchange.com/questions/102933/more-efficient-spatial-join-in-python-without-qgis-arcgis-postgis-etc

Regards, Ahad

[King-David]

Mauricio, please try to keep in mind that some libraries may require installation. Which is A nightmare for the government.

:)

On Tue, Sep 12, 2017 at 9:18 PM Ahad notifications@github.com wrote:

Hi Mauricio There are some libraries in python such as "PySAL" "fiona" can do spatial analysis including spatial join. Writing your own code is not that much hard using such libraries.

This thread at Stackexchange may help you.

https://gis.stackexchange.com/questions/102933/more-efficient-spatial-join-in-python-without-qgis-arcgis-postgis-etc

Regards, Ahad

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/TravelModellingGroup/TMGToolbox/issues/31#issuecomment-329029080, or mute the thread https://github.com/notifications/unsubscribe-auth/AJhKb-bOuNpy3EKxIbwkDVcb5NjKCB9fks5shy1cgaJpZM4POxxm .

[malamillo]

Both 'shapely' and 'fiona' are available within the EMME python libraries installation. Having said that, I'm still looking at the approach suggested by Peter of calling the 'Load Attribute from Polygon' tool. This may be the simplest way to automate the task. Thank you both!

[PeterKucirek]

I should also mention that, as of version 4.3.0, Emme provides a Python API to its Data Tables feature. There is even a spatial_join() function which looks like it will do what you are looking for:

file:///C:/Program%20Files/INRO/Emme/Emme%204/Emme-4.3.3/doc/html/api_doc/emme_api_ref/datatable.html#inro.emme.datatable.Table.spatial_join

Definitely try the tool first, but given the age of the tool, the new API call might be more robust in the long run.

[malamillo]

Agreed, Peter.

Although tagging the links with the screenline id through a spatial join is the easy part.

I still need to get the link cardinal direction (NB, EB, SB, WB) to get the directional flows. I found an old EMME user conference paper online that discusses an approach, however I was wondering if someone has been able to crack down this problem as well.

Source: countgis9.pdf

[PeterKucirek]

Given that the network projection is UTM, with the y-axis representing North, it's a straightforward trigonometric calculation to get the orientation in radians based on the arc from the I-node and J-node.

However, there are two tricks to be aware of:

1. Regions like Peel and Hamilton are roughly 45 degrees oriented, so you might need an "offset" link attribute based on geographical location encoded in the node numbering.
2. Long links with excessively curved shapes (such as a curled highway ramp) will probably need manual handling, if there are any counts on them. Fortunately they can be identified easily by their ratio of shape length to arc length.

[malamillo]

import inro.modeller as m
mm = m.Modeller()
eb=mm.emmebank
scenario = eb.scenario('3100')

load_attribute_from_polygon_package = mm.tool('tmg.network_editing.load_attribute_from_polygon')
load_attribute_from_polygon_package.Scenario = scenario
load_attribute_from_polygon_package.EmmeAttributeIdToLoad = '@test'
load_attribute_from_polygon_package.ShapefilePath = r'E:\Data sources\Mapping\GGHM\screenlinemaps\Cordon_count_program_scr2.shp'
load_attribute_from_polygon_package.ShapeFieldIdToLoad = 'GGHM_id'
load_attribute_from_polygon_package.IntersectionOption = 'INTERSECTS'
load_attribute_from_polygon_package.InitializeAttribute= True
load_attribute_from_polygon_package.run()

---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-15-b338eaf4a269> in <module>()
      8 load_attribute_from_polygon_package.IntersectionOption = 'INTERSECTS'
      9 load_attribute_from_polygon_package.InitializeAttribute= True
---> 10 load_attribute_from_polygon_package.run()
     11 

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-726ed09e989a11e79fae463500000031.pyc in run(self)
    287 
    288         try:
--> 289             self._Execute()
    290 
    291         except Exception, e:

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-726ed09e989a11e79fae463500000031.pyc in _Execute(self)
    303                                      attributes=self._GetAtts()):
    304 
--> 305             polygons = self._LoadPolygons()
    306             print "Loaded polygons"
    307 

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-726ed09e989a11e79fae463500000031.pyc in _LoadPolygons(self)
    372                     polygon[self.ShapefileFieldIdToLoad] = float(val)
    373                 except:
--> 374                     raise Exception("Cannot cast '%s' to float in field '%s' for FID=%s" %(val, self.ShapefileFieldIdToLoad, fid))
    375 
    376                 polygons.append(polygon)

UnboundLocalError: local variable 'val' referenced before assignment

Any idea why this error is being raised?

[PeterKucirek]

Yo dawg, we put an error inside your error so you can err while you err.

Jokes aside, it looks like the field "GGHM_id" does not exist in your shapefile's DBF, which is causing a KeyError than is being mistaken for a TypeError. And then a second error is occurring while the code is trying to handle the error.

[malamillo]

Haha...

good one. Was that an Inception meme?

If "GGHM_id" is a string, would it still be uploaded to '@test'? Would it raise an error?

[PeterKucirek]

Not sure what you're asking. "GGHM_id" is already a string and completely independent in function to '@test'. The former specifies the shapefile attribute to join, while the latter specifies an Emme attribute to save the join into. Input vs. output.

It looks like the field "GGHM_id" does not exist in your shapefile. Python is case-sensitive so double-check that you've got the exact field name.

[malamillo]

The field type for "GGHM_id" in the shapefile is defined as a string, whereas '@test' is defined as a float value.

[PeterKucirek]

That might be OK, as long as the strings are reasonably parsed to numeric values. For example "5.2" is OK, but "five-point-two" is not.

The relevant conversion occurs on L372

[malamillo]

---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-5-df42a663ad0c> in <module>()
      8 load_attribute_from_polygon_package.IntersectionOption = 'INTERSECTS'
      9 load_attribute_from_polygon_package.InitializeAttribute= True
---> 10 load_attribute_from_polygon_package.run()
     11 

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-07fcd0ee98ae11e7b18a463500000031.pyc in run(self)
    287 
    288         try:
--> 289             self._Execute()
    290 
    291         except Exception, e:

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-07fcd0ee98ae11e7b18a463500000031.pyc in _Execute(self)
    303                                      attributes=self._GetAtts()):
    304 
--> 305             polygons = self._LoadPolygons()
    306             print "Loaded polygons"
    307 

C:/Users/ALAMIL~1/AppData/Local/Temp/load_attribute_from_polygon-07fcd0ee98ae11e7b18a463500000031.pyc in _LoadPolygons(self)
    358 
    359     def _LoadPolygons(self):
--> 360         with Shapely2ESRI(self.ShapefilePath) as reader:
    361 
    362             polygons = []

C:/Users/ALAMIL~1/AppData/Local/Temp/geometry-0807f48098ae11e7a349463500000031.pyc in __enter__(self)
    685 
    686     def __enter__(self):
--> 687         self.open()
    688         return self
    689 

C:/Users/ALAMIL~1/AppData/Local/Temp/geometry-0807f48098ae11e7a349463500000031.pyc in open(self)
    438         '''
    439         if self._canread:
--> 440             self._load()
    441         else:
    442             self._create()

C:/Users/ALAMIL~1/AppData/Local/Temp/geometry-0807f48098ae11e7a349463500000031.pyc in _load(self)
    463             type = info[0]
    464             try:
--> 465                 field = self._dbf2fieldMap[type](info[1], info[2], info[3])
    466                 self.fields[field.name] = field
    467             except KeyError, ke:

C:/Users/ALAMIL~1/AppData/Local/Temp/geometry-0807f48098ae11e7a349463500000031.pyc in __init__(self, name, length, decimals, default)
    220             raise IOError("DBF field definition failed: Field length must be at least 3.")
    221         if decimals >= (length - 2):
--> 222             raise IOError("DBF field definition failed: Cannot assign more than {0} decimals for a field length of {1}".format((length - 2), length))
    223 
    224         self.name =str(name)

IOError: DBF field definition failed: Cannot assign more than 11 decimals for a field length of 13

??

[PeterKucirek]

Argg. This error looks like a result of ArcGIS violating their own rules for their shapefile.

At this point I recommend trying it out using Emme's Data Table API. It will be more forgiving.

[malamillo]

Yiaks! Anything else that could be done to address this before jumping into the Data table API stream?

[PeterKucirek]

I'm going to close this issue, but feel free to continue the discussion.