[gdal-dev] densify line?

Bryan Keith bryan at ideotrope.org
Thu Jul 10 14:13:42 EDT 2008

```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

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
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
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
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
```