Numpty's Progress

Tag: awk

  • Convert Manifold aux files to ESRI world files

    The great thing about geospatial metadata standards is that there are so many to choose from. When slinging files about from user to user, it seems that mostly the ESRI stuff (well, the stuff we can pick apart, anyway) is used, plus some sedimentary layers deposited as good ideas from other systems. Commercial GIS systems are much less liberal, because all of them have the inertia of user investment.

    Manifold does things its own way. Rather than using broadly standard prj or wld files, Manifold accompanies every exported file with an XML metafile in its own format. I get sent a lot of bitmaps from Manifold, and if I want to do more than view them on the screen, I need to convert the metadata.

    Below is a small Awk program to convert a Manifold XML metadata file for a bitmap into an ESRI World file:

    #!/bin/awk -f
# convert Manifold xml metadata to ESRI world file
# nb: this is not a general-purpose XML parser
# scruss - 20120508

BEGIN {
    indata = 0;
}

/<coordinateSystem>/ {
    indata = 1;			# in metadata section of file
}

/<\/coordinateSystem>/ {
    indata = 0;			# no longer in metadata
}

/</ && (indata == 1) {
    # build array of items
    sub(/<\/[^>]*>/, "", $0);	# remove end tag
    tagloc=match($0, /<[^>]*>/);
    if ((RSTART == 0) && (RLENGTH == -1)) { 
	# no match!
    }
    else {
	tag=substr($0, RSTART, RLENGTH);
	sub(/</, "", tag);	# remove < and > from tag
	sub(/>/, "", tag);
	value=substr($0, RSTART + RLENGTH); # value is rest of string
	params[tag] = value;
    }
}

END {
    # output six values for ESRI world file
    # all the "+ 0.0" stuff is just to force numeric output
    OFMT="%.10f";
    print params["localScaleX"] + 0.0;
    print params["rotationX"] + 0.0; # these will almost always be zero
    print params["rotationY"] + 0.0;
    print (params["localScaleY"] + 0.0) * -1.0; # negate Manifold value
    print params["localOffsetX"] + 0.0;
    print params["localOffsetY"] + 0.0;
}

    Note that I’m doing what you’re supposed never to do — parse XML as if it were just some text format. If Manifold change their file format even slightly, this program will fail. The above seems to work on all the data I’ve thrown at it, although none of those test files included rotation terms. I don’t think I’ve ever actually seen a world file with rotation defined, so this may be a minor limitation.

  • Drawing Eggs

    When designing wind farms, you need to keep the turbines a certain distance apart from one another. If you don’t, the wakes from the turbines reduce efficiency, and the turbulence can reduce the warranted life of the machine. Typically, a manufacturer might specify a minimum downwind separation of 5 diameters, and a crosswind separation of 3 diameters. It’s an easy check with a buffer overlap, but these buffers are elliptical, which not all GIS packages can draw.

    Take, for example, the following three points designated as (completely made-up)  wind turbine locations:

    name    XCOORD           YCOORD
1    557186.675000    4757125.590000
2    557447.931000    4756968.690000
3    557664.999000    4756817.810000

    These look quite far apart, even if you were using large, 100+m diameter wind turbines:

    But if we have a wind direction of 210°, downwind/crosswind separation of 5D & 3D respectively, and a 101m diameter rotor, it’s not so good:

    Turbines 2 & 3 are too close together; the ellipses shouldn’t touch.

    As awk is the only scripting language I have on my work computer, I wrote the script that generates the buffer shapefile in awk. The script calls Frank’s Shapefile C Library utilities to actually make the shapefile. Here’s the code:

    #!/bin/awk -f
# draw an ellipse based on turbine location to generate
#  for WTG separation buffer
# scruss - 2011-09-27

# assumes that stdin has three columns:
#  1 - label
#  2 - x coordinate
#  3 - y coordinate

# variables:
#  diameter = rotor diameter
#  cross = crosswind separation, diameters
#  down = downwind separation, diameters
#  wind = prevailing wind direction
#  base = base for shape file name

BEGIN {
	OFMT="%.1f";
	CONVFMT="%.1f";
	OFS=" ";

	if (diameter < 0) {
		print "diameter must be set";
		exit;
	}

	if (cross < 0) {
		print "cross must be set";
		exit;
	}

	if (down < 0) {
		print "down must be set";
		exit;
	}

	if (down < cross) {
		print "down must be greater than cross";
		exit;
	}

	if (wind < 0) {
		print "wind must be set";
		exit;
	}

	if (base ~ /^$/) {
		print "base must be a string";
		exit;
	}

	pi = 3.141592654; # I know, I know ...
	# calc cartesian angle from wind bearing, in radians
	beta = ((450 - wind)%360) * pi/180;

	# output shapelib tools init commands
	print "dbfcreate " base " -s name 40";
	print "shpcreate " base " polygon";
}

# for every line
{
	name=$1;
	x=$2;
	y=$3;

	major = diameter * down/2;
	minor = diameter * cross/2;
	first="";
	points="";
	maxn=36;
	for (i=0; i<maxn; i++) {
	    alpha = (i * (360/maxn)) * pi/180;
	    x1 = x + major * cos(alpha) * cos(beta) - minor * sin(alpha) * sin(beta);
	    y1 = y +  major * cos(alpha) * sin(beta) + minor * sin(alpha) * cos(beta);
	    if (i == 0) { # store the first point
		first= x1 " " y1;
	    }
	    points = points  " " x1 " " y1;
	}
	points = points  " " first;
	print "dbfadd " base ".dbf " name;
	print "shpadd " base, points;
}

    awk is charmingly odd in that you can specify variable on the command line. Here’s how I called it, with the above coordinates as input:

    awk -v diameter=101 -v cross=3 -v down=5 -v wind=210 -v base="fakewtg-ellipse" -f separation.awk

    Pipe the output through a shell, and there are your ellipses.