def line2pts(geom, dl=None):
    """Given an input line geom, generate points at fixed interval
    
    Useful for extracting profile data from raster
    """
    #Extract list of (x,y) tuples at nodes
    nodes = geom.GetPoints()
    #print "%i nodes" % len(nodes)
   
    #Point spacing in map units
    if dl is None:
        nsteps=1000
        dl = geom.Length()/nsteps

    #This only works for equidistant projection!
    #l = np.arange(0, geom.Length(), dl)

    #Initialize empty lists
    l = []
    mX = []
    mY = []

    #Add first point to output lists
    l += [0]
    x = nodes[0][0]
    y = nodes[0][1]
    mX += [x]
    mY += [y]

    #Remainder
    rem_l = 0
    #Previous length (initially 0)
    last_l = l[-1]
    
    #Loop through each line segment in the feature
    for i in range(0,len(nodes)-1):
        x1, y1 = nodes[i]
        x2, y2 = nodes[i+1]
      
        #Total length of segment
        tl = np.sqrt((x2-x1)**2 + (y2-y1)**2)

        #Number of dl steps we can fit in this segment
        #This returns floor 
        steps = int((tl+rem_l)/dl)

        if steps > 0:
            dx = ((x2-x1)/tl)*dl
            dy = ((y2-y1)/tl)*dl
            rem_x = rem_l*(dx/dl)
            rem_y = rem_l*(dy/dl)
            
            #Loop through each step and append to lists
            for n in range(1, steps+1):
                l += [last_l + (dl*n)]
                #Remove the existing remainder
                x = x1 + (dx*n) - rem_x
                y = y1 + (dy*n) - rem_y
                mX += [x]
                mY += [y]

            #Note: could just build up arrays of pX, pY for entire line, then do single z extraction
            #Update the remainder
            rem_l += tl - (steps * dl)
            last_l = l[-1]
        else:
            rem_l += tl 

    return l, mX, mY
Given an input line geom, generate points at fixed interval
    
    Useful for extracting profile data from raster