For a layer's feature to be represented in all layers it must be represented in any selection of two layers a well. Using this theory you can simply run multiple intersection functions according to the order below. The example assumes 4 total layers named Layer1, Layer2, Layer3, and Layer4.
- Run an intersection between Layer1 and Layer2 to get Output1.
- Run an intersection between Output1 and Layer3 to get Output2.
- Run an intersection between Output2 and Layer4 to get Output3.
Output 3 should contain all features that appear in all 4 layers. For the 'Count' field, I think you can just file it with the number of layers involved in the intersect, unless I am mistaking what data you want that field to hold.
This method isn't eloquent and if you have many layers or layers with many features it could prove time consuming but it should work and is a simple process to put together.
Your script is too complex, simply use the modulo function:
def azimuth(point1, point2):
'''azimuth between 2 shapely points (interval 0 - 360)'''
angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
return np.degrees(angle) if angle >= 0 else np.degrees(angle) + 360
azimuth(interP,P2)
112.61986494834154
azimuth(P2,interP)
292.61986494834156
azimuth(P1,P2)
207.64597536482526
azimuth(P2,P1)
27.645975364825276
azimuth(P1,P2)
207.64597536482526
azimuth(P2,P1)
27.645975364825276
# result
azimuth(P2,interP) - azimuth(P1,interP)
12.550802249443507
# or
azimuth(interP,P2) - azimuth(interP,P1)
12.550802249443507
But you use here the Euclidean distance and if your points are geodetic (angles, longitude, latitude) this distance is meaningless (see inaccurate distance measurements in Python).
With the pygc module (Vincenty's formulae)
print P1.distance(P2) # shapely distance in "angular unit"
2.37065391828006e-05
from pygc import * # distance in meters
print great_distance(start_latitude=P1.y, start_longitude=P1.x, end_latitude=P2.y,end_longitude=P2.x)
{'reverse_azimuth': 27.553194547658386, 'distance': array(2.619663014755096), 'azimuth': 207.55319613451832}
P2az = great_distance(start_latitude=P2.y, start_longitude=P2.x, end_latitude=interP.y,end_longitude=interP.x)['reverse_azimuth']
P1az = great_distance(start_latitude=P1.y, start_longitude=P1.x, end_latitude=interP.y,end_longitude=interP.x)['reverse_azimuth']
print P2az -P1az
12.592170938638404
And as Charles Petzold concludes (How Far from True North are the Avenues of Manhattan? ), the angular difference is minimal
Best Answer
Instead of Intersect, use Clip.
Vector > Geoprocessing Tools > Clip