The Problem We have to filter out the roads and ditches without removing streams that cross roads or follow them closely. I’m going to use PostGIS to find the intersection of the streams lines data with a buffered roads polygon. If the intersected line is less than 50% of the length of the stream line, […]
Posts Tagged ‘Database’
Posted by smathermather on June 17, 2016
Posted by smathermather on May 29, 2016
Here was the problem that needed solved last week (we have a few similar problems in upcoming projects, so this was an exciting thing to try): we needed to use PostGIS to access data in a SQLServer database. The SQLServer database backs the web site in question, the underlying content management system, etc., so no– removing SQLServer isn’t really an option at this stage. Obviously PostGIS is a requirement too… .
Before I go further, I used tds_fdw as the foreign data wrapper. The limitations here are as follows: it is a read-only connection, and only works with Sybase and SQLServer, as it uses tabular data stream protocol for communicating between client and server. This is not as generic a solution as we can use. Next time I’ll try ogr_fdw which is more generic as it can connect with other databases and other data types. Another advantage to ogr_fdw is we can use IMPORT FOREIGN SCHEMA. Regina Obe and Leo Hsu warn though to limit this to 150 tables or so for performance reasons.
With the limitations listed above, this is how we built the thang:
DROP SERVER beach_fdw CASCADE; -- Create the server connection with FOREIGN DATA WRAPPER CREATE SERVER beach_fdw FOREIGN DATA WRAPPER tds_fdw OPTIONS (servername 'name_or_ip', port '1433', database 'non-spatial-database', tds_version '7.1', msg_handler 'notice'); -- We map the postgres user to the user that can read the table we're interested in CREATE USER MAPPING FOR postgres SERVER beach_fdw OPTIONS (username 'user', password 'password'); -- Create the actual foreign table connection CREATE FOREIGN TABLE beach_closures ( AutoNumber int NOT NULL, Title varchar NOT NULL, StartDate timestamp without time zone NOT NULL, WaterQuality varchar NOT NULL, Latitude varchar NOT NULL, Longitude varchar NOT NULL, BeachStatus varchar NOT NULL, ClosureColor varchar NOT NULL) SERVER beach_fdw OPTIONS (schema_name 'schema_name', table_name 'vw_CMPBeachClosures'); -- Now we create a spatial view using our longitude and latitude CREATE VIEW v_beach_closures AS SELECT AutoNumber, Title, StartDate, WaterQuality, Latitude, Longitude, BeachStatus, ClosureColor, ST_SetSRID(ST_MakePoint(Longitude::numeric, Latitude::numeric), 4326) AS geom FROM beach_closures;
Voila! A nice little PostGIS enabled view of a SQLServer view or table!
Posted by smathermather on January 4, 2016
Trying out PostGIS 2.2’s ST_ApproximateMedialAxis capabilities today. I’m going to use it to label parcels. So here is my parcel fabric:
And here’s what the skeleton of that fabric looks like:
It get’s pretty interesting where the parcels are more complicated:
Oh, the code you say? Well, it couldn’t be much easier:
SELECT gid, ST_ApproximateMedialAxis(geom) AS geom FROM cuy_parcel_2015;
Or the full version for a new table:
DROP TABLE IF EXISTS cuy_parcel_2015_medial; CREATE TABLE cuy_parcel_2015_medial AS ( SELECT gid, ST_ApproximateMedialAxis(geom) AS geom FROM cuy_parcel_2015 );
ALTER TABLE cuy_parcel_2015_medial ADD PRIMARY KEY (gid); CREATE INDEX cuy_parcel_2015_medial_geom_idx ON cuy_parcel_2015_medial USING GIST ( geom ); Happy New Year.
Posted by smathermather on November 15, 2014
This is a quick blog post about technologies that I don’t know well… so please comment if you know better. GeoGig and dat are great tools for addressing versioning in data, so what’s the difference?
GeoGig is built on Java and meant for any “simple features” geometry (points, lines, polygons).
It’s strength is that it is built from the ground up to handle geometries well, going beyond CRUD functions to specifically address geospatial problems in versioning. Think of it as git for geospatial data.
From the website:
“Users are able to import raw geospatial data (currently from Shapefiles, PostGIS or SpatiaLite) in to (sic) a repository where every change to the data is tracked. These changes can be viewed in a history, reverted to older versions, branched in to sandboxed areas, merged back in, and pushed to remote repositories.”
Ok, so how about dat?
“Dat is an open source project that provides a streaming interface between every file format and data storage backend.”
A cursory look indicates it will work for geospatial data, but effectively as blobs, with no special handling for changes within features like GeoGig. But, it does what GeoGig does not, and that is to make datasets automatically syncable.
Like all projects, each has its strengths. Choose your project wisely.
Posted by smathermather on January 9, 2014
Free book to the first to tell me what the is name of the children’s book on the armchair in the foreground:
Book editing in progress… .
Posted by smathermather on August 6, 2012
Anyone who has been following along at home knows I don’t think much like a DBA. Sometimes that’s good; mostly it’s probably bad. In this post, I hope it will be interesting.
The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM. This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions. In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention. This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust. Doing one well seems to preclude doing the other well.
Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts. Here’s a decent list, which I haven’t checked for timeliness: http://vterrain.org/Elevation/contour.html. Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now. I tried GRASS, but it was too slow for the amount of data being processed. The command in ArcGIS ends up looking something like this:
TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #
And thus we go from this:
Now to do it 623 more times… . And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.
Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:
Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:
So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles. An ST_Touches should do the trick here:
SELECT p.tile, a.tile AS tiles FROM cuy_contour_tiles AS p, cuy_contour_tiles As a WHERE ST_Touches(p.geom, a.geom);
But, that does not arrange our table very nicely:
I want something that looks like this:
N2225655, N2225660, N2220650 …
N2225660, N2225665, N220665 … etc.
Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.
Adapting their entry, we can reorganize those entries like this:
SELECT STRING_AGG(a.tile, ',' ORDER BY a.tile) As tiles FROM cuy_contour_tiles AS p, cuy_contour_tiles As a WHERE ST_Touches(p.geom, a.geom) GROUP BY p.tile ORDER BY p.tile;
And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:
SELECT p.tile, 'TopoToRaster_sa "' || 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;' || 'W:\contours\cuy_contours\shapefiles\' || STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile) || '.shp elevation Contour;" ' || '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2' || ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"' || ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # ' As tiles FROM cuy_contour_tiles AS p, cuy_contour_tiles As a WHERE ST_Touches(p.geom, a.geom) GROUP BY p.tile ORDER BY p.tile
Ok. I have 623 more tiles to do. Until next time… .
Posted by smathermather on May 17, 2012
First we create a function to do our angle calculations, then use select to loop through all the records and do the calculations.
Within our function, first we find our first (in this case) five nearest streets using our distance operator , based on bounding boxes, then we further order by ST_Distance, with a LIMIT 1 to get just the closest street. By the way, this is about twice as fast as the last approach.
And our new code:
CREATE OR REPLACE FUNCTION angle_to_street (geometry) RETURNS double precision AS $$ WITH index_query as (SELECT ST_Distance($1,road.geom) as dist, degrees(ST_Azimuth($1, ST_ClosestPoint(road.geom, $1))) as azimuth FROM street_centerlines As road ORDER BY $1 <#> road.geom limit 5) SELECT azimuth FROM index_query ORDER BY dist LIMIT 1; $$ LANGUAGE SQL; DROP TABLE IF EXISTS address_points_rot; CREATE TABLE address_points_rot AS SELECT addr.*, angle_to_street(addr.geom) FROM address_points addr;
I have to credit Alexandre Neto with much of the code for this solution.
Posted by smathermather on May 12, 2012
My favorite cartography book is Eduard Imhof’s Cartographic Relief Presentation. A few years back I picked this book up (translated to English) from ESRI press for $75 if memory serves me. Now it can be gotten for much cheaper.
Imhof spends a lot of time on feature simplification and separation, a problem which keeps me up at night. For example, if you have a lot of overlapping and/or touching contours, what are the local distortions and simplifications that can be done to enhance map readability?
This problem applies to other overlapping features, such as a road following a river, where at a given scale one might not distinguish between the two, so we exaggerate their differences, in this case, while maintaining their basic topology.
I’ve played with this problem before, but now I’m considering a new approach. I pinged Martin Davis of JTS fame for some advice. He suggested using force directed layout to solve this problem, much as he does in JTS Test Builder for magnifying/investigating topology.
Posted by smathermather on May 10, 2012
In a previous post (I feel like I say that a lot), I wrote about rotating address points to match nearby roads in replicate the effect of USGS quads that represented small buildings with little squares that followed the nearby road alignment.
The function was effective:
ALTER TABLE jstein.address_points_all ADD COLUMN azimuth integer; UPDATE jstein.address_points_all addr SET azimuth = degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom))) FROM jstein.street_centerlines road;
but deadly slow when applied to all 500,000 address points. And so we iterate. First, I’ll show you our reasonably fast, but elegantly written solution. Later I’ll have a post on our very fast, but somewhat hackerish approach.
A reminder of our target aesthetic:
And our new code:
CREATE TABLE address_points_rot AS WITH index_query AS ( SELECT addr.*, ST_Distance(addr.the_geom, road.the_geom) AS distance, degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom))) as azimuth FROM address_pointssub as addr, street_centerlines as road WHERE ST_DWithin(addr.the_geom, road.the_geom, 2000) = true ORDER BY addr.the_geom <-> road.the_geom ) SELECT DISTINCT ON(cartodb_id) * FROM index_query ORDER BY gid, distance
In the interest of full disclosure, I stole the structure of this from Paul Ramsey’s post on K Nearest Neighbor searches in PostGIS. AFAIK, you need Postgresql 9.1 and Postgis 2.0 or later to do this… . For 500,000 points this took 16 minutes. Narrowing that search distance would help a lot. I suspect there is a much better way to structure this using true KNN for each point, but my SQL-fu failed me here… . Suggestions welcome.
Posted by smathermather on May 10, 2012
For one of our applications, we need 8-point compass roses placed at each of our points, as well as a circle 40 meters in diameter as well as one 140 meters in diameter. We did a bit of work with this a while back in GeoServer using SLDs. Now we’d like to refine it, and implementing this in an SLD is beyond my skills.
So, we move to the back end. The following is a little function to construct the roses for us, which we’ll then display though GeoServer. Just an FYI, this is running on a PostGIS 1.3 database, and I didn’t (at the time anyway) grok the use of ST_RotateZ, so there are some explicit position calls, rather than just creating the crosshairs and rotating it:
CREATE OR REPLACE FUNCTION pie ( geom Geometry, filling numeric, crust numeric ) RETURNS Geometry AS $$ DECLARE slices Geometry; BEGIN --slices = ST_MakeLine(geom, ST_Translate(geom, dough, pan)); slices = ST_Union( ST_Collect( ST_Collect( ST_Collect( ST_MakeLine(ST_Translate(geom, filling, 0), ST_Translate(geom, crust, 0)) ,ST_MakeLine(ST_Translate(geom, 0, filling), ST_Translate(geom, 0, pi()/4 * crust)) ) ,ST_Collect( ST_MakeLine(ST_Translate(geom, -filling, 0), ST_Translate(geom, pi()/4 * -crust, 0)) ,ST_MakeLine(ST_Translate(geom, 0, -filling), ST_Translate(geom, 0, pi()/4 * -crust)) ) ) ,ST_Collect( ST_Collect( ST_MakeLine(ST_Translate(geom, filling * pi()/4, filling * pi()/4), ST_Translate(geom, crust * pi()/4, crust * pi()/4)) ,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, filling * pi()/4), ST_Translate(geom, -crust * pi()/4, crust * pi()/4)) ) ,ST_Collect( ST_MakeLine(ST_Translate(geom, filling * pi()/4, -filling * pi()/4), ST_Translate(geom, crust * pi()/4, -crust * pi()/4)) ,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, -filling * pi()/4), ST_Translate(geom, -crust * pi()/4, -crust * pi()/4)) ) ) ) ,ST_Collect( ST_ExteriorRing(ST_Buffer(geom, filling, 32)) ,ST_ExteriorRing(ST_Buffer(geom, crust, 32)) ) ) ; RETURN slices; END; $$ LANGUAGE plpgsql;
So now to test it with a single point:
DROP TABLE IF EXISTS test_pie; CREATE TABLE test_pie AS SELECT pie(ST_MakePoint(0,0), 40, 140);
And view in uDig, or application of choice: