1

I have to create a 1.5 km buffer on each latitude/longitude of layer A. I would like to know how many points of line B will be inside the buffer of layer A. I used the following code but it takes too much time.

Is there a faster way to do this?

       crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic
       crs_bng = proj.Proj(init='epsg:4326') 

       # then cast your geographic coordinate pair to the projected system


       d1=OrderedDict()
       d2=OrderedDict()
       temp=[]
       temp1=[]
       lat=0.0
       lon=0.0
       a1=0.0

       templist=[]

       Alatlong=[]
       attrs=A.getFeatures()
       hlidx = A.fieldNameIndex('Y')
       #print idx
       hloidx = A.fieldNameIndex('X')


       Alatlong=  [(i[hlidx],i[hloidx]) for i in attrs]       

       Blatlong=[]
       attrs1=B.getFeatures()
       plidx = B.fieldNameIndex('Latitude')
       #print idx
       ploidx = B.fieldNameIndex('Longitude')

       Blatlong=  [(j[plidx],j[ploidx]) for j in attrs1]
       distance =1500

       for i in Alatlong:
            ans=2

            input_lat=i[0]
            input_lon=i[1]
            x, y = proj.transform(crs_wgs, crs_bng, input_lon, input_lat)

            point_1 = geometry.Point(x, y)
            circle_buffer = point_1.buffer(distance)                               

            for j in Blatlong:                 

                input_lat1=j[0]
                input_lon1= j[1]                 

                x1, y1 = proj.transform(crs_wgs, crs_bng, input_lon1, input_lat1)

                point_2=geometry.Point(x1,y1)                

                x=[]
                if point_2.within(circle_buffer):
                   print 'True'

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.