Archive for the ‘Well Known Text’ Category
Terminal one-liner to create shapefile from WKT using ogr2ogr
Fair warning, this is a Linux-themed solution. But I expect it could be ported to Windows without too much trouble.
Today I needed to create a shapefile with a single point geometry in it to use as an input for a random GIS utility executable, which required its input to be a an actual shapefile rather than an array of coordinate pairs, or even a single coordinate pair, or even separate x=/y= input parameters for a single coordinate pair. I already had the point coordinate I wanted to use, and I didn’t want to go to the trouble to make this shapefile! So I got to wondering if I could use ogr2ogr to render a simple shapefile from the terminal.
Fortunately, with the help of some utilities built in to my Ubuntu shell, I was able to come up with the following solution. This compound command 1) creates an empty file (dataset.csv), 2) pushes two features into it, and 3) uses ogr2ogr to translate the CSV into a shapefile projected accordingly–in this case, to EPSG:4326.
Here’s the full command. Note that a pair of ampersands (&&) are used to daisy-chain the individual steps together into a single instruction:
touch dataset.csv && printf "gid,WKT,some_value\n1,POINT(-82.048051 33.567181),Alpha\n2,POINT(-92.7774 35.9829),Bravo\n" > dataset.csv && ogr2ogr -f "ESRI Shapefile" dataset.shp -dialect sqlite -sql "SELECT gid, GeomFromText(WKT), some_value FROM dataset" dataset.csv -a_srs EPSG:4326
Here’s a breakdown of what happens..
touch dataset.csv
- Create a file (in the current directory) named “dataset.csv”
printf "gid,WKT,some_value\n1,POINT(-82.048051 33.567181),Alpha\n2,POINT(-92.7774 35.9829),Bravo\n" > dataset.csv
- Push a string through standard output (STDOUT) and INTO (>) the new file “dataset.csv”
- The string imposes line breaks using the “\n” sequence and establishes a CSV type containing three lines
- The three lines of content include: a column header, feature 1, and feature 2
ogr2ogr -f "ESRI Shapefile" dataset.shp -dialect sqlite -sql "SELECT gid, GeomFromText(WKT), some_value FROM dataset" dataset.csv -a_srs EPSG:4326
- The ogr2ogr command exports a shapefile, “dataset.shp”
- The input file, “dataset.csv” is included toward the end, along with the declaration of the source data’s CRS (EPSG:4326/WGS84)
- ogr2ogr is instructed to use the “SQLITE” dialect of sql to pull data from the target dataset
- The SQL command makes sure ogr2ogr knows to identify the WKT column as the source data’s geometry field
It runs pretty quickly, and does exactly what I wanted: Quickly crank out a shapefile if all I have is a point or a polygon or something, already in WKT format. I’m sure I’ll be using this again, so I wanted to document it here for easy reference later. Hopefully someone else finds it useful!
Convert Google Maps Polygon (API V3) to Well Known Text (WKT) Geometry Expression
There’s dozens of reasons why you might want the Well Known Text (WKT) geometry expression for a Google Maps Polygon object.
Assuming you’re using the Google Maps API V3, and you’ve got a variable referencing your Polygon, I’ll suggest two approaches you can take to iterate over the paths and vertices in your Google Maps polygon and return the geometry expression as a Well Known Text string.
Add a Simple Utility Method to Your Project
Easy enough. Just add the following method to your project. Look below the method for an example of how you’d call it.
function GMapPolygonToWKT(poly)
{
// Start the Polygon Well Known Text (WKT) expression
var wkt = "POLYGON(";
var paths = poly.getPaths();
for(var i=0; i<paths.getLength(); i++)
{
var path = paths.getAt(i);
// Open a ring grouping in the Polygon Well Known Text
wkt += "(";
for(var j=0; j<path.getLength(); j++)
{
// add each vertice and anticipate another vertice (trailing comma)
wkt += path.getAt(j).lng().toString() +" "+ path.getAt(j).lat().toString() +",";
}
// Google's approach assumes the closing point is the same as the opening
// point for any given ring, so we have to refer back to the initial point
// and append it to the end of our polygon wkt, properly closing it.
//
// Also close the ring grouping and anticipate another ring (trailing comma)
wkt += path.getAt(0).lng().toString() + " " + path.getAt(0).lat().toString() + "),";
}
// resolve the last trailing "," and close the Polygon
wkt = wkt.substring(0, wkt.length - 1) + ")";
return wkt;
}
Here’s how you’d access the Well Known Text expression using the utility method:
// Assuming you've already instantiated "myPolygon" somewhere.
var wkt = GMapPolygonToWKT(myPolygon);
Extend Google’s Polygon Object Prototype with a ToWKT() Method
There’s nothing wrong with the first approach, but you might find it handy to extend Google’s Polygon object prototype, itself, to include a ToWKT() member function, which makes it even easier to get its Well Known Text. To do that, add the following JavaScript somewhere near the top of your code (caveat—this will need to be called after the Google Maps library has been loaded):
if (typeof google.maps.Polygon.prototype.ToWKT !== 'function') { google.maps.Polygon.prototype.ToWKT = function() { var poly = this; // Start the Polygon Well Known Text (WKT) expression var wkt = "POLYGON("; var paths = poly.getPaths(); for(var i=0; i<paths.getLength(); i++) { var path = paths.getAt(i); // Open a ring grouping in the Polygon Well Known Text wkt += "("; for(var j=0; j<path.getLength(); j++) { // add each vertice, automatically anticipating another vertice (trailing comma) wkt += path.getAt(j).lng().toString() + " " + path.getAt(j).lat().toString() + ","; } // Google's approach assumes the closing point is the same as the opening // point for any given ring, so we have to refer back to the initial point // and append it to the end of our polygon wkt, properly closing it. // // Additionally, close the ring grouping and anticipate another ring (trailing comma) wkt += path.getAt(0).lng().toString() + " " + path.getAt(0).lat().toString() + "),"; } // resolve the last trailing "," and close the Polygon wkt = wkt.substring(0, wkt.length - 1) + ")"; return wkt; }; }
If you prefer the second approach, you can get the Well Known Text expression like this:
// Assuming you've already instantiated "myPolygon" somewhere.
var wkt = myPolygon.ToWKT();
PostGIS: query all multipolygon parcels with at least one hole
I was writing some code to iterate over Well Known Text expressions for polygon features, and I decided I needed to test the most complex edge-case I could think of–multipolygon geometries where at least one of the bound polygons has a hole (i.e. an interior ring).
I ended up with the following query. This seems like the kind of thing I’ll want to reuse later, so I’m noting it here. For good measure, I also use a rudimentary technique to sort the output with the most complicated geometries in the table at the top of the list. Basically, the more “text” it takes to describe the geometry using Well Known Text, the larger and more complex I figure it must be!
SELECT SomePrimaryId, /* your primary key, i.e. ogc_fid, etc. */ SomeUniqueId, /* your descriptive id, i.e. a parcel number */ ST_NumGeometries(wkb_geometry) AS num_geoms, ST_NRings(wkb_geometry) AS num_rings, ST_AsText(ST_Centroid(wkb_geometry)) AS center, Char_Length(ST_AsText(wkb_geometry)) AS len, ST_AsText(wkb_geometry) AS wkt FROM SomePolygonTable WHERE ST_NumGeometries(wkb_geometry) > 1 AND ST_NRings(wkb_geometry) > ST_NumGeometries(wkb_geometry) ORDER BY Char_Length(ST_AsText(wkb_geometry)) ASC ;
Just for the sake of promoting caution, I’m not certain this is a definitive approach for identifying the largest geometry in a table, as the length of the binary representation and the length of the readable text representation do not correspond one-to-one. Moreover, a feature could have more vertices that required less precision to express (fewer decimal position), than a geometry with fewer vertices that needed more precision, and then you have to ask, which is bigger, fewer vertices and more text, or more vertices that coincidentally did not require as much text? My conclusion is, the “most complicated geometry” is probably relative to the one asking the question. However for my purposes, this was close enough to put the most complicated stuff at the top of the list.