[gdal-dev] densify line?

[Bryan Keith]

Matt,

I did something like this a while ago in ogr and Python.  Some of the code
is attached.  You may find it useful.  It's certainly not in the form of a
ready-to-use utility.

Bryan

> Hello All,
>
> Does ogr have a densify line option or utility?
>
> The goal is to densify (add vertices at regular intervals) to a 1
> dergree by 1 degree fishnet grid before projecting from WGS84 to UTM (or
> whatever) so the result follows earth curvature rather then jumping
> straight from corner point to corner point.
>
> I could use Arcinfo DENSIFYARC, but would like to avoid having to
> convert to a coverage if I can.
>
> thanks,
>
> --
> matt wilkie
> --------------------------------------------
> Geographic Information,
> Information Management and Technology,
> Yukon Department of Environment
> 10 Burns Road * Whitehorse, Yukon * Y1A 4Y9
> 867-667-8133 Tel * 867-393-7003 Fax
> http://environmentyukon.gov.yk.ca/geomatics/
> --------------------------------------------
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
import ogr
import math
#AddVerticesToLN will add vertices to a LN (mLN) so that no vertices will be
#farther than dDist from each other
#instead of creating the LN, this routine will just return a list of vertices
#currently this routine only supports LNs where mLN.GetGeometryCount() == 0
def AddVerticesToLN(mLN, dDist):
  lPTs = []
  for i in range(mLN.GetPointCount() - 1):
    #get the current two vertices and their distance
    X1 = mLN.GetX(i)
    Y1 = mLN.GetY(i)
    X2 = mLN.GetX(i + 1)
    Y2 = mLN.GetY(i + 1)
    dDXY = ReturnDistance(X1, X2, Y1, Y2)
    #include the current vertex in the output list
    lPTs += [[X1, Y1]]
    #if the distance between the vertices is greater than dDist, then I
    #need to add vertices
    if (dDXY > dDist):
      #how many segments do I need to make?
      fCount = float(dDXY) / dDist
      iCount = int(fCount)
      if (iCount != fCount):
        iCount += 1
      #calculate the distance of each new segment
      dDistNew = dDXY / iCount
      dDistAlong = dDistNew
      #make the new PTs
      for j in range(iCount - 1):
        lPTs += [AlongSimple(X1, Y1, X2, Y2, dDistAlong)]
        dDistAlong += dDistNew
  #add the final vertex
  lPTs += [[X2, Y2]]
  return lPTs
#AlongSimple is a simple version of the along function since I don't have a
#current need to implement a real version.
#This simply takes 2 XY pairs as input and a distance along the LN that
#they create.
#It returns the new XY at dDist along the LN
def AlongSimple(X1, Y1, X2, Y2, dDist):
  #vertical LN case
  if (X1 == X2):
    X3 = X1
    if (Y1 < Y2):
      Y3 = Y1 + dDist
    else:
      Y3 = Y1 - dDist
  #horizontal LN case
  elif (Y1 == Y2):
    Y3 = Y1
    if (X1 < X2):
      X3 = X1 + dDist
    else:
      X3 = X1 - dDist
  #diagonal LN case
  else:
    m = (Y2 - Y1) / (X2 - X1)
    a = math.atan(m)
    dX = math.cos(a) * dDist
    dY = math.sin(a) * dDist
    if (X1 < X2):
      Y3 = Y1 + dY
      X3 = X1 + dX
    else:
      Y3 = Y1 - dY
      X3 = X1 - dX
  return [X3, Y3]
def ReturnDistance(x1, x2, y1, y2, z1 = 0, z2 = 0):
  return ((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2) ** 0.5