Welcome to Planet PostGIS

November 27, 2014

Boston GIS (Regina Obe, Leo Hsu)

osm2pgrouting for Windows

For loading OSM data in a format already ready for pgRouting queries, the two common tools I've seen used are osm2po and osm2pgrouting. Unfortunately osm2pgrouting has been for a long time a Unix only tool mostly because no one tried to compile it for windows or test it on windows to see if it works. So this means that windows users who wanted a quick road to pgRouting nirvana had to use osm2po. Turns out the osm2pgrouting code compiles fine on windows (at least under mingw-w64 chain) and seems to work fine, so we should start making it available for windows.

This is still experimental, but if you are a windows user (or have a pgRouting windows user friend) and want to kick the tires a bit, you'll find osm2pgrouting for windows (both 32-bit and 64-bit) binaries on the PostGIS windows download page in the Unreleased section.

I compiled osm2pgrouting (for pgRouting 2.0) against PostgreSQL 9.4 using the same mingw-w64 chain I use to build PostGIS and pgRouting for windows -- clumsily detailed in my gist shell scripts - . Though it depends on libpq, it seems to work fine if you use it against the PostgreSQL 9.3 libpq. So just need to copy the additional files (osm2pgrouting and expat) into your PostgreSQL bin folder and you are good to go loading data into a pgRouting enabled database. It is assumed you already have PostGIS 2.1+ and pgRouting 2.0 already installed which come bundled together as part of the PostGIS 2.1+ windows installers and also available on the PostGIS windows download page.

Which should you use: osm2pgrouting or osm2po?

In many ways osm2po is superior to osm2pgrouting and has had more time invested in it. It can handle larger files on smaller machines, can work with the OSM pbf format, and it supports both Unix like and windows platforms. That said osm2po is not absolutely superior to osm2pgrouting. osm2po is inferior in a couple of ways that are important to me.

  • It's not open source, it's freeware -- which means you can't kick the code around and learn from it or inspect it. It also means you can't change it and contribute back your changes.
  • The primary motive of osm2po is not as a loading tool for pgRouting. osm2po is a complete routing engine with webserver which doesn't even rely on PostgreSQL. That said, if the maintainer decides one day providing a pgRouting export is no longer beneficial to osm2po development path, or pgRouting devs need to tweak it a bit to work with newer pgRouting, it will be difficult.
  • It's written in Java. This is my own personal dislike. Not that anything is wrong with Java, it's just sometimes a bit of a nuisance to get the VM going on a server I happen to get stuck on that I can't install anything I want. osm2pgrouting is much lighter weight not requiring a JVM, though it does require PostgreSQL,PostGIS, and pgRouting -- which you should have already.
  • osm2po generates an intermediary sql file which you then have to load via psql. I prefer the direct load approach of osm2pgrouting over having a big-o sql file to contend with after.

That said osm2po is great in many ways, but alternatives are always nice to have. Only annoying thing with having alternatives is deciding which alternative is best for what you are doing.

by Regina Obe (nospam@example.com) at November 27, 2014 11:34 PM

November 21, 2014

Postgres OnLine Journal (Leo Hsu, Regina Obe)

PostGIS Day synopsis

Yesterday was PostGIS day or as some may call it, Post GIS day and a couple of interesting things happened this day:

  • PostgreSQL 9.4 RC1 came out.
  • There were parties and unconferences, many summarized on http://2014.postgisday.rocks
  • I managed to entertain myself with a Conway's game of life PostGIS raster Map algebra style and pondered how wonderful PostGIS would be if it could generate animated gifs with some sort of aggregate function; to which I was politely called crazy by some of my fellow PSC friends.

But what was greatest of all and took the cake were these pictures:


Continue reading "PostGIS Day synopsis"

by Leo Hsu and Regina Obe (nospam@example.com) at November 21, 2014 09:16 PM

November 20, 2014

Boundless Geo

Happy PostGIS Day!

PostGISYesterday was GIS Day, which means that today is PostGIS Day — get it? Post-GIS Day! To celebrate, we’re giving a 50% discount on online PostGIS training through the end of the week! Visit our training catalog and use promo code “postgis2014″ to take advantage of this offer.

A lot has happened since last year, when I joined Stephen Mather, Josh Berkus, and Bborie Park in a PostGIS all-star hangout with James Fee.

In case you missed them, here are some features from our blog and elsewhere that highlight what’s possible with PostGIS:

An, as always, be sure to check out our workshops for a slew of PostGIS-related courses, including Introduction to PostGIS and our Spatial Database Tips & Tricks.

Interested in learning about or deploying PostGIS? Boundless provides support, training, and maintenance for installations of PostGIS. Contact us to learn more.

The post Happy PostGIS Day! appeared first on Boundless.

by Paul Ramsey at November 20, 2014 01:54 PM

Boston GIS (Regina Obe, Leo Hsu)

PostGIS Day Game of Life celebration

This year's PostGIS day, I decided to celebrate with a little Conway's Game of Life fun inspired by Anita Graser's recent blog series Experiments with Game of Life. The path I chose to simulate the Game of life is a little different from Anita's. This variant exercises PostGIS 2.1+ raster mapalgebra and PostgreSQL 9.1+ recursive queries. Although you can do this with PostGIS 2.0, the map algebra syntax I am using is only supported in PostGIS 2.1+. My main disappointment is that because PostGIS does not yet support direct generation of animated gifs I had to settle for a comic strip I built unioning frames of rasters instead of motion picture. Hopefully some day my PostGIS animated gif dream will come true.


Continue reading "PostGIS Day Game of Life celebration"

by Regina Obe (nospam@example.com) at November 20, 2014 07:04 AM

November 18, 2014

Postgres OnLine Journal (Leo Hsu, Regina Obe)

FOSS4GNA 2015 and PGDay San Francisco March 2015

PGDay 2015 San Francisco will be held March 10th 2015 in Hyatt, San Franciso Airport, Burlingame, CA (Just outside of San Francisco). This year PGDay will be hosted along-side Free and Open Source Geospatial North America (FOSS4GNA) conference 2015 which runs March 9th-12th and EclipseCon NA 2015. Speaker submissions for FOSS4GNA 2015 and EclipseCon NA 2015 will end this Monday November 17th, 2014.


Continue reading "FOSS4GNA 2015 and PGDay San Francisco March 2015"

by Leo Hsu and Regina Obe (nospam@example.com) at November 18, 2014 03:46 AM

November 03, 2014

Boston GIS (Regina Obe, Leo Hsu)

AvidGeo Conference: LocationTech Tour Boston, Dec 8 2014

We'll be giving a PostGIS Introduction tutorial at AvidGeo Conference: LocationTech Tour Boston which will be held 8 December 2014 from 8:30 AM to 4:30 PM at:

Hack / Reduce
275 Third St
Cambridge, MA, 02142

Tickets are $40 for admission.

In addition to our PostGIS Intro tutorial, there is a great line-up of talks from visionaries such as Adena Schutzberg, talking: State of the Geospatial Industry and Raj Singh on the topic of The NoSQL Geo Landscape. Also check out other tutorials on CartoDb and QGIS.

by Regina Obe (nospam@example.com) at November 03, 2014 06:40 AM

October 26, 2014

Stephen Mather

Airspace — A deep rabbit hole

In previous maps we looked at Class B, C, and D airspace. Let’s add in Class E0 and E5… (not yet in 3D):

(Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL)

Map showing Class B, C, D, E0, and E5 airspace

Map showing Class B, C, D, E0, and E5 airspace

Previous posts:

http://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again/

and

http://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again-reprise/


by smathermather at October 26, 2014 02:16 AM

October 25, 2014

Stephen Mather

Airspace is complicated — and so I abuse PostGIS once again — Reprise…

In the previous post: http://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again/ we explore the 3D shape and complexity of controlled airspace.

Now here’s the rest of the code. We’ll add our affine transformation ala Seth Fitsimmons:

    SELECT 
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )
    AS geom

And integrate that into our original function:

-- Inputs are a geometry and an elevation to move the geometry to
CREATE OR REPLACE FUNCTION threed_iso(footprint geometry, elevation numeric)
  RETURNS geometry AS
$BODY$

-- Force 3D, then translate to the input elevation
WITH floor AS
(
    SELECT ST_Translate( ST_Force3DZ(footprint), 0, 0, elevation ) AS geom
),
-- Now make isometric (and begin Seth Code)
iso AS
(
    SELECT 
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )

    AS geom
    FROM floor
)
-- We'll force it back to 3D so QGIS is happy
SELECT ST_Force2D(geom) FROM iso
;
$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

And voila!

DROP TABLE IF EXISTS class_c_isoc;

CREATE TABLE class_c_isoc AS
	SELECT gid, airspace, name, lowalt, highalt, threed_iso(geom, lowalt::numeric * 5) AS geom
	FROM class_c_subset;

Let’s take a look at Washington, DC and surrounds, another nice complicated example:

3D Figure of DC controlled airspace

3D Figure of DC controlled airspace

And again with map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL:

3D Figure of DC controlled airspace

3D Figure of DC controlled airspace

 


by smathermather at October 25, 2014 05:17 PM

Stephen Mather

Airspace is complicated — and so I abuse PostGIS once again

Let’s ignore for a moment the drone hobbiest / enthusiast. What is the shape of airspace for airplanes and commercial and government unmanned aircraft flying under Certificates of Authorization, and how can we visualize it?

Thanks to Anita in the last post, we have the Class B,C,D,E Airspace Shape Files which helps us define the overall shape of controlled airspace.

Map of Detroit, Cleveland, and Pittsburgh Class B Airspace

Map of Detroit, Cleveland, and Pittsburgh Class B Airspace

But, these are 3D things. I want to visualize them thus. Let us put some constraints on the problem. Let’s do it all in PostGIS, that way we can see it in QGIS or on the web. Let’s not use full PostGIS 3D (i.e. the SFCGAL library), not because it isn’t awesome (it truly is) but because it can be hard to install at the moment (although see https://github.com/vpicavet/docker-pggis for an easy way with docker). Finally, true 3D with unconstrained viewing angles and 100% flexibility is… is… well it usually sucks. So, we’ll stick to isometric viewing (thanks to Seth Fitzsimmons of Stamen http://stamen.com/ for his PostGIS isometric code which will be released upon his word). (Update — all the code is there…):

-- Inputs are a geometry and an elevation to move the geometry to
CREATE OR REPLACE FUNCTION threed_iso(footprint geometry, elevation numeric)
  RETURNS geometry AS
$BODY$

-- Force 3D, then translate to the input elevation
WITH floor AS
(
    SELECT ST_Translate( ST_Force3DZ(footprint), 0, 0, elevation ) AS geom
),
-- Now make isometric (and begin Seth Code)
iso AS
(
    SELECT
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )

    AS geom
    FROM floor
)
-- We'll force it back to 3D so QGIS is happy
SELECT ST_Force2D(geom) FROM iso
;
$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

Ok, now let’s rock some geometries with this bad function:

DROP TABLE IF EXISTS class_c_isoc;

CREATE TABLE class_c_isoc AS
	SELECT gid, airspace, name, lowalt, highalt, threed_iso(geom, lowalt::numeric * 5) AS geom
	FROM class_c_subset;

And what do our controlled airspaces look like?

Isometric view of Cleveland controlled airspace

Isometric view of Cleveland controlled airspace

Kind of conical, in this case with some “wings” that serve as approaches. It makes sense, I guess. At the center, where the airport is, controlled airspace goes from the ground up. As we get further from the airport, a set of concentric rings starting at higher elevations provide a controlled space that allows for landing and taking off.

Image of Cleveland airspace with Stamen Toner map as backdrop

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

There are more complicated ones, naturally. We need look no further than Detroit for additional complexity:

Visualization of Detroit controlled airspace

Visualization of Detroit controlled airspace

airspace_detroit_toner

No wonder every time I fly through Detroit, no matter where I’m going, I end up flying over Lake Erie.

If we want really complicated, we need look no further than Cincinnati:
airspace_cinci

What is going on there? Is that due to shape of the hills, or city boundaries?

Finally, what does airspace look like over Ohio, West Virginia, and Pennsylvania (etc.), overall?
airspace_is_complicated_regional

And while the following map isn’t quite right, here is a figure including many of the small airports sans controlled airspace:

View of all controlled and uncontrolled airspace across Ohio and neighbors.

View of all controlled and uncontrolled airspace across Ohio and neighbors.

May an aeronautical pox be upon you!
The above line was not intended in bad taste, but just an homage to the “red dot fever” of early neo-geography (neo-geography which I’m informed definitionally doesn’t exist). Only a few minutes ago, it dawned on me that the deleted phrase could be misinterpreted these days… .

On a related note, if you want an interesting analysis of ebola, read Zeynep Tufekci’s analysis.

(for the record, all heights are exagerated by 5x, for clarity of reading).

Also, in case it wasn’t obvious: Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.


by smathermather at October 25, 2014 06:13 AM

October 20, 2014

Oslandia

Oslandia releases "Cuardo", 3D GIS viewer


Oslandia announces today the preliminary result of one of its research and development project : Cuardo .



Cuardo is an OpenSource WebGL 3D data viewer, focusing on urban data analysis and visualization.



Cuardo features :

  • 3D terrain data model
  • Multi-textured terrain (satellite, thematic)
  • Draping of 2D vector data
  • Textured 3D meshes
  • Custom styling, with 3D effects (extrusion…)
  • Level of detail management
  • Classic 3D user navigation
  • WebServices integration ( WMS, WFS )
  • Identify tool with custom actions (popup, links…)
  • Focus on performances
  • High customization possibilities

The specificity of this tool is that it deals really well with 3D features as objects, including LOD management, whereas most 3D viewers cannot distinguish objects and make it hard for object picking. Being able to visualize a full city in your browser without any plugin to install, and interact with other client application, makes Cuardo a good option for application development requiring 3D visualization, which you can mix with other great web libraries.


Cuardo is on GitHub, get it, fork it, PR it : https://github.com/Oslandia/cuardo


This software component is part of our 3D software stack. This Cuardo release is a first public release, and it is still under heavy development, but already usable for interesting use cases, as you can see in the video above.


We actively work on a full-featured 3D software solution to focus on 3D geographical data management. Our aim is to provide a full set of components allowing to :

  • Import 3D data
  • transform 2D data into 3D data
  • Store high volumes of data
  • Perform 3D analysis
  • Combine 3D data with - elevation models - point clouds - Raster data - 2D vector data - alphanumeric data
  • Setup webservices for data diffusion, including 3D data
  • Visualize the data in a browser

The release of Cuardo web client is the missing piece for a vertical 3D GIS data management stack.


Oslandia has a strong expertise on 3D GIS, be it on the server side, or as stated by this anouncement, on the client side.


Oslandia is commited to improving this software solution, adapting it and making it fit your specific use cases. Do not hesitate to get in touch if you want some more information on how to take advantage of your 3D spatial data in the best way. We will be delighted to hear your use cases : infos+3d@oslandia.com .


We would like to thank the Grand Lyon metropolitan area, for distributing a lot of GIS data as opendata, including textured 3D models : http://www.smartdata.grandlyon.com

by Vincent Picavet at October 20, 2014 08:00 AM

October 09, 2014

Boston GIS (Regina Obe, Leo Hsu)

FOSS4G NA 2015 Call for Speakers now Open

FOSS4G NA 2015 (Free and Open Source Geospatial North America conference) will be running March 9th-12th 2015 In Burlingame, California (Hyatt, San Francisco Airport). Call for talks and workshops is now open and we've already got some submissions to review. Deadline for submission is November 17th 2014. Presenters including (workshop and those giving 35-minute standard talk) get to attend the 3-days general sessions for free. The conference consists of one day of workshops and three days of general sessions, with lightning talks on some or all of of the general session days.

Details about presentations here: CFP details. The FOSS4G NA 2015 is a collaborative event co-sponsored by OSGeo and LocationTech

FOSS4G NA 2015 will be co-hosted with EclipseCon North America 2015. Admission to FOSS4G NA 2015 gives you admission to EclipseCon NA 2015. So get 2 conferences for the price of one and go conference hopping.

For PostgreSQL and PostGIS folks there is a likely chance that there will be a PostgreSQL day at the same venue during the same period. So we might have a very special PostGIS / PostgreSQL day. I'll update as more details of that surface.

by Regina Obe (nospam@example.com) at October 09, 2014 03:02 AM

October 01, 2014

Boundless Geo

PostGIS Training: Creating Overlays

PostGIS training At Boundless, we recognize that software is only as good as the problems it can be used to solve. That’s why we don’t just develop the best open source spatial software, we also teach you how to use it with supporttraining, and workshops directly from the experts.

One question that comes up often during our PostGIS training is “how do I do an overlay?” The terminology can vary: sometimes they call the operation a “union” sometimes an “intersect”. What they mean is, “can you turn a collection of overlapping polygons into a collection of non-overlapping polygons that retain information about the overlapping polygons that formed them?”

screenshot_67

So an overlapping set of three circles becomes a non-overlapping set of 7 polygons. screenshot_69

Calculating the overlapping parts of a pair of shapes is easy, using the ST_Intersection() function in PostGIS, but that only works for pairs, and doesn’t capture the areas that have no overlaps at all. How can we handle multiple overlaps and get out a polygon set that covers 100% of the area of the input sets? By taking the polygon geometry apart into lines, and then building new polygons back up. Let’s construct a synthetic example: first, generate a collection of random points, using a Gaussian distribution, so there’s more overlap in the middle. The crazy math in the SQL below just converts the uniform random numbers from the random() function into normally distributed numbers.

CREATE TABLE pts AS
WITH rands AS (
  SELECT generate_series as id, random() AS u1, random() AS u2 FROM generate_series(1,100)
)
SELECT
  id,
  ST_SetSRID(ST_MakePoint(
    50 * sqrt(-2 * ln(u1)) * cos(2*pi()*u2),
    50 * sqrt(-2 * ln(u1)) * sin(2*pi()*u2)),4326) AS geom
FROM rands;

The result looks like this: screenshot_70 Now, we turn the points into circles, big enough to have overlaps.

CREATE TABLE circles AS
SELECT id, ST_Buffer(geom, 10) AS geom FROM pts;

Which looks like this: screenshot_71 Now it’s time to take the polygons apart. In this case we’ll take the exterior ring of the circles, using ST_ExteriorRing(). If we were dealing with complex polygons with holes, we’d have to use ST_DumpRings(). Once we have the rings, we want to make sure that everywhere rings cross the lines are broken, so that no lines cross, they only touch at their end points. We do that with the ST_Union() function.

CREATE TABLE boundaries AS
SELECT ST_Union(ST_ExteriorRing(geom)) AS geom
FROM circles;

What comes out is just lines, but with end points at every crossing. screenshot_72 Now that we have noded lines, we can collect them into a multi-linestring and feed them to ST_Polygonize() to generate polygons. The polygons come out as one big multi-polygon, so we’ll use ST_Dump() to convert it into a table with one row per polygon.

CREATE SEQUENCE polyseq;
CREATE TABLE polys AS
SELECT nextval('polyseq') AS id, (ST_Dump(ST_Polygonize(geom))).geom AS geom
FROM boundaries;

Now we have a set of polygons with no overlaps, only one polygon per area. screenshot_73 So, how do we figure out how many overlaps contributed to each incoming polygon? We can join the centroids of the new small polygons with the set of original circles, and calculate how many circles contain each centroid point. screenshot_74 A spatial join will allow us to calculate the number of overlaps.

ALTER TABLE polys ADD COLUMN count INTEGER DEFAULT 0;
UPDATE POLYS set count = p.count
FROM (
  SELECT count(*) AS count, p.id AS id  
  FROM polys p 
  JOIN circles c 
  ON ST_Contains(c.geom, ST_PointOnSurface(p.geom)) 
  GROUP BY p.id
) AS p
WHERE p.id = polys.id;

That’s it! Now we have a single coverage of the area, where each polygon knows how much overlap contributed to it. Ironically, when visualized using the coverage count as a variable in the color ramp, it looks a lot like the original image, which was created with a simple transparency effect. However, the point here is that we’ve created new data, in the count attribute of the new polygon layer. screenshot_75 The same decompose-and-rebuild-and-join-centroids trick can be used to overlay all kinds of features, and to carry over attributes from the original input data, achieving the classic “GIS overlay” workflow. Happy geometry mashing!

Want to learn more? Try our Introduction to PostGIS online training course!

The post PostGIS Training: Creating Overlays appeared first on Boundless.

by Paul Ramsey at October 01, 2014 03:33 PM

September 30, 2014

Boston GIS (Regina Obe, Leo Hsu)

Waiting for PostGIS 2.2 - ST_ClipByBox2D - Map dicing Redux in concert with ST_Tile

One of the new features coming in PostGIS 2.2 is ST_ClipByBox2D (thanks to Sandro Santilli's recent commits funded by CartoDB ). However to take advantage of it, you are going to need your PostGIS compiled with GEOS 3.5+ (very recent build) which has not been released yet. Windows folks, the PostGIS 2.2 9.3 and 9.4 experimental binaries are built with the latest GEOS 3.5 development branch, so you should be able to test this out with Winnie's experimental builds.

Since the dawn of PostGIS, PostGIS users have needed to mutilate their geometries in often painful and horrific ways. Why is ST_ClipByBox2D function useful, because its a much faster way of mutilating your geometries by a rectangular mold than using ST_Intersection. Why would you want to mutilate your geometries? There are many reasons, but I'll give you one: As your geometry approaches the area of your bounding box and as your bounding box size decreases, the more efficient your spatial index becomes. You can consider this article, Map dicing redux of the article I wrote (eons ago) - Map Dicing and other stuff which describes the same approach with much older technology. Though I will be using more or less the same dataset Massachusetts TOWNSSURVEY_POLYM (its newer so you can't really compare) and I tried to simplify my exercise a bit (not resorting to temp tables and such), my focus in this article will be to compare the speed between the new ST_ClipByBox2D approach and the old ST_Intersection approach. The spoiler for those who don't have the patience for this exercise is that using ST_ClipByBox2D at least on my sample data set on my puny Windows 7 64-bit desktop using PostgreSQL 9.4beta2 64-bit was about 4-5 times faster than using ST_Intersection. There was a downside to this speedup. With the ST_Intersection approach, I had no invalid polygons. In the case of ST_ClipByBox2D, I had one invalid polygon. So as noted in the docs, use with caution. We'd be very interested in hearing other people's experiences with it.

One other benefit that ST_ClipByBox2D has over ST_Intersection which I didn't test out is that although ST_Intersection doesn't work with invalid geometries, ST_ClipByBox2D can.

For these exercises, I'm going to also abuse the PostGIS raster function ST_Tile for geometry use, cause for some strange reason, we have no ST_Tile function for PostGIS geometry. For those who missed the improper use of ST_Tile for raster, refer to Waiting for PostGIS 2.1 - ST_Tile.


Continue reading "Waiting for PostGIS 2.2 - ST_ClipByBox2D - Map dicing Redux in concert with ST_Tile"

by Regina Obe (nospam@example.com) at September 30, 2014 02:55 PM

September 26, 2014

PostGIS Development

Getting distinct pixel values and pixel value counts of a raster

PostGIS raster has so so many functions and probably at least 10 ways of doing something some much much slower than others. Suppose you have a raster, or you have a raster area of interest — say elevation raster for example, and you want to know the distinct pixel values in the area. The temptation is to reach for ST_Value function in raster, but there is a much much more efficient function to use, and that is the ST_ValueCount function.

ST_ValueCount function is one of many statistical raster functions available with PostGIS 2.0+. It is a set returning function that returns 2 values for each row: a pixel value (value), and a count of pixels (count) in the raster that have that value. It also has variants that allow you to filter for certain pixel values.

This tip was prompted by the question on stackexchange How can I extract all distinct values from a PostGIS Raster?

Continue Reading by clicking title hyperlink ..

by Regina Obe at September 26, 2014 12:00 AM

September 25, 2014

Archaeogeek (Jo Cook)

New Portable GIS Releases

I’m pleased to announce not one, but two new releases of Portable GIS! The first, version 4.2, contains QGIS 2.4 and PostGIS 1.5 and will be the last release to include that version of PostGIS.

The second, version 5, contains QGIS 2.4 and PostGIS 2.1 and all future releases will be based on this. Get them here:

There are two important things about these two releases. The first is that the download size is almost twice what it was previously- this is out of my control and is due to the increased size of the QGIS installation. The second is that version 5.0 should be considered a beta release as it’s very much un-tested. If there are any problems let me know via the Portable GIS google group.

You might also be interested to see a short presentation I did about Portable GIS at GeoMob London in September: http://archaeogeek.github.io/geomob-portablegis/#/

September 25, 2014 08:49 AM

September 24, 2014

Oslandia

Back from FOSS4G 2014 - part 2



Apart from Oslandia's participation at FOSS4G 2014 in Portland ( See part 1), a lot of topics have been discussed during this week of conference.


Presentation's slides and video are being put online and you can replay the whole conference if you could not attend, but here is a really short summary of a few things we noticed during the five days.

  • Fiona and Rasterio are now mature modern GIS Python libraries on top of GDAL/OGR, and are continuously improved.
  • WebGL is all around the place. With new support from the various browser's editors, this is definitly the backend of choice for large data rendering on online maps. This is still the exploration phase though, and there is currently no definite solution for webgl mapping, as things are evolving really fast and use cases appear.
  • Use cases for Opensource GIS software continue to grow larger. It even goes to space with PlanetLabs satellite constellation !
  • OpenLayers 3 and Leaflet are the two client libraries of choice, with distinct use cases, but an excellent quality for both
  • Design is here, for the good sake of our eyes and user experience. This is the first time that a specific track is fully dedicated to maps design.
  • Vector strikes back ! With the evolution of browser's technologies, vector data is gaining momentum. Vector tiles are a natural evolution for this, and formats and standards start to emerge. More performance, less space, and a door open for new end-user features
  • Following the vector movement, symbology also evolves to become dynamic. What if those colors nicely blend from one to another, and the thickness of these lines are controlled by a function of the zoom level ? Better UX !
  • GIS involves large dataset, and opensource software are now ready to crunch these teràbytes of data. Be it for mosaïcking rasters for the whole world, machine learning to predict climate change, or distributed data analysis, data size limits continue to be pushed further.
  • Well-known OpenSource components continuously gain features : Mapserver, QGIS, GDAL/OGR, PostGIS improve from version to version
  • CartoDB raised 8 million USD and continues to grow
  • A large part of the Mapbox team was present and showed various technologies which are at the bleeding edge of online mapping
  • A lot of various large organization were present, and we can notice TriMet, US DoD, Amazon (sponsor), GitHub, ESRI, RedHat… Big players are in the place too

FOSS4G is technical and at the same time user-oriented. It gathers the opensource community and industrials as well as academics. This is the best way to foster an intense brainstorming and create new ideas, visions and technologies.


Do not hesitate to check the videos, and give us your favorites, twitting @oslandia or mail us on any of these subjects at infos@oslandia.com !

by Vincent Picavet at September 24, 2014 05:00 PM

September 23, 2014

Postgres OnLine Journal (Leo Hsu, Regina Obe)

FOSS4G and PGOpen 2014 presentations

At FOSS4G we gave two presentations. The videos from other presentations are FOSS4G 2014 album. I have to commend the organizers for gathering such a rich collection of videos. Most of the videos turned out very well and are also downloadable in MP4 in various resolutions. It really was a treat being able to watch all I missed. I think there are still some videos that will be uploaded soon. As mentioned lots of PostGIS/PostgreSQL talks (or at least people using PostGIS/PostgreSQL in GIS).


Continue reading "FOSS4G and PGOpen 2014 presentations"

by Leo Hsu and Regina Obe (nospam@example.com) at September 23, 2014 06:55 PM

Oslandia

Back from FOSS4G 2014 - part 1



FOSS4G , the main international OSGeo event, takes place once per year and gather all folks involved in opensource geospatial technologies at large.

This year, the conference took place in Portland, Oregon, and gathered around 1000 people, on very various topics ranging from deeply technical Mapserver pro tips, to burnout management, through opendata policies and map designs.

This year's edition was really good, as always for this conference. A lot of social events allowed direct talks to passionate people from all over the world, and the presentations were diverse and of high quality.

Oslandia at FOSS4G

Oslandia is a full player in the OpenSource GIS field, and we participated this year with various interventions.


First of all, we gave a workshop on our 3D stack, with new features and a now vertical stack of software, from storage to visualization.

This workshop features the use of this stack for managing 3D GIS data. Given an opendata dataset from the city of Lyon, France , we integrate the data into PostGIS, and use PostGIS 3D features to analyze and enrich this data. Then, Horao ( http://oslandia.github.io/horao/ ), a 3D destktop viewer coupled with QGIS, is used to do some 3D visualization.


In a second part of the workshop, we use the very recently released Cuardo web client ( http://www.cuardo.org ) to display the data in a standard browser, using WebGL.



The workshop is self-contained and you can follow it completely to discover these new exciting features in OpenSource 3D GIS. find it online : https://github.com/Oslandia/workshop-3d

Gimme some YeSQL !

Vincent Picavet also did a presentation "Gimme some YeSQL - and a GIS -", which focused on latest and upcoming PostgreSQL features. From updateable views to the new JSONB data type, this is a practical and technical presentation of the big improvements coming in PostgreSQL 9.4. The talk illustrates the use of these features in a GIS context.

Slides are online As well as the video (the title is wrong).



OpenSource software for surveying

Another presentation by Vincent was an insight on a project we achieved for GIS Apavil in Romania, and concerns a full GIS for water distribution network and wastewater networks. The presentation especially shows how we setup some specific opensource surveying tools, and the way we deal with offline, versioning and synchronization of data.

Slides are on GitHub and the video should be up soon.

3D GIS

Olivier Courtin showed a summary of the current state of the 3D stack we develop and its foreseen future. Slides and video are online, and there should be more information coming soon.


Slides : https://github.com/Oslandia/presentations/raw/master/foss4g_2014/gisgoes3d_olivier_courtin.pdf



Get in touch !

Oslandia will continue to make research and development on all of these topics. If you are interested in any of these subjects, do not hesitate to get in touch at infos@oslandia.com


A post to follow will give you some background information on what we heard at FOSS4G and where to hear more !

by Vincent Picavet at September 23, 2014 05:00 PM

September 22, 2014

Paul Ramsey

PostGIS Feature Frenzy

A specially extended feature frenzy for FOSS4G 2014 in Portland. Usually I only frenzy for 25 minutes at a time, but they gave me an hour long session!

PostGIS Feature Frenzy — Paul Ramsey from FOSS4G on Vimeo.

Thanks to the organizers for giving me the big room and big slot!

by Paul Ramsey (noreply@blogger.com) at September 22, 2014 11:30 PM

September 16, 2014

Paul Ramsey

PostGIS for Managers

At FOSS4G this year, I wanted to take a run at the decision process around open source with particular reference to the decision to adopt PostGIS: what do managers need to know before they can get comfortable with the idea of making the move.

The Manager's Guide to PostGIS — Paul Ramsey from FOSS4G on Vimeo.

by Paul Ramsey (noreply@blogger.com) at September 16, 2014 07:48 PM

September 11, 2014

Postgres OnLine Journal (Leo Hsu, Regina Obe)

FOSS4G 2014 televised live

If you weren't able to make it to FOSS4G 2014 this year, you can still experience the event Live. All the tracks are being televised live and its pretty good reception. https://2014.foss4g.org/live/. Lots of GIS users using PostGIS and PostgreSQL. People seem to love Node.JS too.

After hearing enough about Node.JS from all these people, and this guy (Bill Dollins), I decided to try this out for myself.

I created a node.js web application - which you can download from here: https://github.com/robe2/node_postgis_express . It's really a spin-off from my other viewers, but more raw. I borrowed the same ideas as Bill, but instead of having a native node Postgres driver, I went for the pure javascript one so its easier to install on all platforms. I also experimented with using base-64 encoding to embed raster output directly into the browser so I don't have to have that silly img src path reference thing to contend with.

by Leo Hsu and Regina Obe (nospam@example.com) at September 11, 2014 08:37 PM

September 10, 2014

A GeoSpatial World

Transformation de formats 3D (OBJ, PLY, STL, VTP) vers PostGIS WKT POLYHEDRALSURFACE


3D (OBJ,PLY,STL,VTP) vers WKT POLYHEDRALSURFACE

image

image_thumb175_thumb ntroduction

L’idée de l’outil que je vais vous présenter est venu de la lecture d’un message dans la liste PostGIS-User polyhedral surface import from common 3D formats , l’auteur du message recherche un outil permettant de traduire quelques formats 3D (OBJ, PLY, STL ..) vers le format WKB/WKT pour une géométrie de type POLYHDERAL SURFACE et je me suis dis alors pourquoi ne pas développer un tel outil, le viewer 3D déjà réalisé me permettant de valider la transformation.
J’ai développé cet outil en C# 2013 avec :
  • ActiViz .Net un wrapper pour VTK
  ActiViz .Net


vtk-logo_thumb2 VTK
  • ScintillaNET un puissant contrôle d’édition de textes pour Windows et un wrapper pour le composant Scintilla.
image ScintillaNET
Scintilla
J’ai ensuite recherché sur internet des fichiers exemple pour chaque format de données 3D que j’ai inclus dans le répertoire Data de l’installation de l’outil. Voici quelques liens ou vous trouverez des fichiers à télécharger :
  • OBJ
OBJ Files
  • PLY
PLY Files
  • STL
STL Files Sharing
 

image image
image image image

image_thumb174_thumb nstallation


image_thumb41 Cliquez sur le lien pg3DImport_setup.exe pour télécharger le programme d’installation, puis lancer l’exécution.
Cliquer sur le  bouton Suivant à chaque étape, puis sur le bouton Terminer.
image
image
image
image
image
image

image_thumb179_thumb tilisation


image Double cliquez sur l’icône créé par le programme d’installation pour lancer l’application.
image
image

Ouvrir et transformer un fichier 3D

image Cliquez sur cet icône pour choisir un fichier 3D à transformer.
Allez dans le répertoire d’installation de Pg3DImport puis le sous-répertoire DATA, puis sélectionnez le type de fichier 3D que vous souhaitez transformer OBJ Files par exemple puis sélectionnez le fichier cesna.obj et cliquez sur le bouton Ouvrir.

image
Vous pourrez alors visualiser le fichier dans la partie viewer.

image
image Cliquez sur cet icône pour obtenir la transformation vers le format POLYHEDRALSURFACE WKT.


image

Vous pouvez alors copier le résultat dans la partie éditeur pour l’insérer dans une table PostGIS (version 2.0 ou supérieur) ou bien visualiser le résultat dans l’application Pg3DViewer que vous avez certainement déjà installée.

image
image
Visualisation dans Pg3DViewer

image

Un autre exemple

Ouverture d’un fichier STL dans le répertoire DATA de l’installation.
image
Affichage dans Pg3DViewer.

image
Vous pouvez ouvrir et convertir les différents fichiers présents dans le répertoire DATA de l’installation, mais si vous avez des fichiers 3D, utilisez plutôt ceux-ci.

image_thumb100_thumb1 onclusion.

Nous sommes arrivés à la fin de ce billet, un pas de plus dans le nouveau monde 3D de PostGIS, n’hésitez pas à me faire part de vos idées d’amélioration.

by Jérôme ROLLAND (noreply@blogger.com) at September 10, 2014 09:38 AM

A GeoSpatial World

Transformation de formats 3D (OBJ, PLY, STL, VTP) vers PostGIS WKT POLYHEDRALSURFACE

 

3D (OBJ,PLY,STL,VTP) vers WKT POLYHEDRALSURFACE

image

image_thumb175_thumb ntroduction

L’idée de l’outil que je vais vous présenter est venu de la lecture d’un message dans la liste PostGIS-User polyhedral surface import from common 3D formats , l’auteur du message recherche un outil permettant de traduire quelques formats 3D (OBJ, PLY, STL ..) vers le format WKB/WKT pour une géométrie de type POLYHDERAL SURFACE et je me suis dis alors pourquoi ne pas développer un tel outil, le viewer 3D déjà réalisé me permettant de valider la transformation.

J’ai développé cet outil en C# 2013 avec :

  • ActiViz .Net un wrapper pour VTK
  ActiViz .Net


vtk-logo_thumb2 VTK
  • ScintillaNET un puissant contrôle d’édition de textes pour Windows et un wrapper pour le composant Scintilla.
image ScintillaNET
  Scintilla

J’ai ensuite recherché sur internet des fichiers exemple pour chaque format de données 3D que j’ai inclus dans le répertoire Data de l’installation de l’outil. Voici quelques liens ou vous trouverez des fichiers à télécharger :

  • OBJ

OBJ Files

  • PLY

PLY Files

  • STL
STL Files Sharing

 

image image  
image image image

 

image_thumb174_thumb nstallation

 

image_thumb41

Cliquez sur le lien pg3DImport_setup.exe pour télécharger le programme d’installation, puis lancer l’exécution.

Cliquer sur le  bouton Suivant à chaque étape, puis sur le bouton Terminer.

  image
  image
  image
  image
  image
  image

 

image_thumb179_thumb tilisation

 

image Double cliquez sur l’icône créé par le programme d’installation pour lancer l’application.
  image
  image

Ouvrir et transformer un fichier 3D

image Cliquez sur cet icône pour choisir un fichier 3D à transformer.
  Allez dans le répertoire d’installation de Pg3DImport puis le sous-répertoire DATA, puis sélectionnez le type de fichier 3D que vous souhaitez transformer OBJ Files par exemple puis sélectionnez le fichier cesna.obj et cliquez sur le bouton Ouvrir.

image
  Vous pourrez alors visualiser le fichier dans la partie viewer.

image
image Cliquez sur cet icône pour obtenir la transformation vers le format POLYHEDRALSURFACE WKT.


image

Vous pouvez alors copier le résultat dans la partie éditeur pour l’insérer dans une table PostGIS (version 2.0 ou supérieur) ou bien visualiser le résultat dans l’application Pg3DViewer que vous avez certainement déjà installée.

 
image
  image
  Visualisation dans Pg3DViewer

image

Un autre exemple

  Ouverture d’un fichier STL dans le répertoire DATA de l’installation.
image
  Affichage dans Pg3DViewer.

image

Vous pouvez ouvrir et convertir les différents fichiers présents dans le répertoire DATA de l’installation, mais si vous avez des fichiers 3D, utilisez plutôt ceux-ci.

 

image_thumb100_thumb1 onclusion.

Nous sommes arrivés à la fin de ce billet, un pas de plus dans le nouveau monde 3D de PostGIS, n’hésitez pas à me faire part de vos idées d’amélioration.

by Jérôme ROLLAND (noreply@blogger.com) at September 10, 2014 07:55 AM

PostGIS Development

PostGIS 2.1.4 Released

The 2.1.4 release of PostGIS is now available.

The PostGIS development team is happy to release patch for PostGIS 2.1, the 2.1.4 release. As befits a patch release, the focus is on bugs, breakages, and performance issues

http://download.osgeo.org/postgis/source/postgis-2.1.4.tar.gz

Continue Reading by clicking title hyperlink ..

by Regina Obe at September 10, 2014 12:00 AM

September 09, 2014

Safe Software

Free and Open Source: The Inside Scoop on FME and FOSS

“Open source tools play well with others.” Paul Ramsey, PostGIS founder, talks Open Source at the FME International User Conference 2014.

“Open source tools play well with others.” Paul Ramsey, PostGIS founder, talks Open Source at the FME International User Conference 2014.

As a software developer, coding something that’s already been designed, implemented, and optimized is a pretty gigantic waste of time.

Open source libraries help developers avoid re-inventing wheels. At this year’s FME UC, PostGIS founder Paul Ramsey took the stage to talk about the awesomeness that is Open Source. He explained how it’s important that developers use available knowledge—and FME devs understand this, adding value by building useful functionality on top of open source libraries.

Safe Software proudly supports the Free and Open Source Software (FOSS) community and regularly sponsors the North American FOSS4G conferences. FME uses many FOSS libraries, and when the devs here at Safe aren’t playing soccer or LEGO or a game of oversized Jenga, they’re working hard to make these libraries more accessible.

Our use of these libraries is a friendly two-way street. A two-way Yellow Brick Road, if you will. FOSS adds value to FME’s advanced data transformation and automation capabilities, while our contributions to these technologies have included funding, coding, testing, and bug fixes. Here are ten such libraries.

10 FOSS Libraries that FME Uses and Contributes to

This Mario-like map is made possible by the integration of Mapnik and FME.

This Mario-like map is made possible by the integration of Mapnik and FME.

1. Mapnik

Mapnik is a toolkit that renders raster and vector data into beautiful cartographic maps. FME’s Mapnik transformer lets users leverage Mapnik’s capabilities with any of the 325+ data formats FME supports. With Mapnik and FME, you can also create Web Map Tiling Services, make a map out of your face, map yer way to buried treasure arr matey!, and other such highly critical tasks. FME’s intuitive graphical interface makes it easy for any user to benefit from this complex library.

2. SpatiaLite

SpatiaLite extends SQLite into a simple Spatial DBMS. FME’s SpatiaLite Reader / Writer allows users to connect this format with hundreds of other formats, and leverage all of FME’s data transformation and automation capabilities. For example, you can instantly import SpatiaLite data into Google Maps Engine using our free Data Loader for Google Maps Engine.

FME reads and writes a variety of open source databases.

FME reads and writes a variety of open source databases.

3. PostGIS

PostGIS is a spatial database that adds support for geographic objects and location queries to PostgreSQL. FME’s PostGIS Reader / Writer lets users integrate hundreds of data sources with this format. Watch Paul Ramsey’s FME UC talk to see the evolution of PostGIS, from the ideas that inspired it to its widespread use today.

4. JTS and GEOS

JTS and GEOS (Geometry Engine – Open Source) provide spatial predicates and functions for processing 2D geometry. JTS is the Java library and GEOS is the C++ library. Some of FME’s powerful geometry transformations are based on GEOS/JTS algorithms.

5. OGR and GDAL (Geospatial Data Abstraction Library)

These OSGeo libraries support reading/writing over 50 raster formats (GDAL) and 20 vector formats (OGR). FME uses these libraries to add several readers and writers to its list of supported formats, which means these formats gain a wide range of enterprise transformation and automation capabilities.

Bonus: if you use a GDAL-based format in FME, you can feel comfortable knowing it was probably implemented by our in-house world Tetris champion.

6. FDO

FDO works with a variety of geospatial data. This OSGeo API brings additional formats to FME, for example SQLite. In addition, Safe offers the free FME FDO Provider to allow AutoCAD Map 3D users to read 9 additional data formats. Safe was one of the earliest adopters of this API (version 3.0) and is acknowledged as an official Product Using FDO.

7. Ingres

Ingres is an RDBMS that supports large commercial and government applications. Safe’s partnership with Actian resulted in an open-source FME reader/writer for Ingres tables. This means increased data accessibility for Ingres users, making it easy to integrate spatial and enterprise data.

FME Workbench and the FME Data Inspector are available on Linux.

FME Workbench and the FME Data Inspector are available on Linux.

8. Qt

Qt is an application framework that we’ve used to help make FME fully cross-platform. Since FME 2013, FME has been available in beta on Mac and Linux platforms.

9. Interlis

Interlis is a geodata modeling and integration language adopted by the Swiss government as their national standard. Eisenhut Informatik is one of Safe’s longest standing open source partners, and their ili2fme plug-in allows Interlis data to be translated to/from any of the formats supported by FME.

10. LASzip

LASzip is a specialized library by Martin Isenburg that compresses LiDAR data without information loss—turning chubby LAS point clouds into lean, manageable LAZ files. The incorporation of this technology into FME has resulted a great improvement in the ability to archive and exchange point cloud data.

But Wait, There’s More

For the full list of libraries licensed under the LGPL, visit our FOSS Code Center. You can also view more information about FME’s open source libraries by opening FME Desktop and selecting About > Legal Notices.

FME for Free

Data Loader for Google Maps Engine

The Data Loader for Google Maps Engine is one free tool we offer for converting data with FME technology.

FME is FOSS! Sort of. We do offer some free data translation tools.

In addition to the FME FDO Provider mentioned above, we also offer the Data Loader for Google Maps Engine, the CityGML Importer for InfraWorks, and the FME Revit Exporter Plug-in.

The Easy Translator is also available as a free web service, for immediately translating data into your required format and coordinate system.

Moving forward, we aim to make more free tools and web services like these.

banner_fmerocks

How have you integrated your software or project with Free and Open Source Software? Are you attending FOSS4G this year, or have you in the past? Let us know in the comments!

The post Free and Open Source: The Inside Scoop on FME and FOSS appeared first on Safe Software Blog.

by Tiana Warner at September 09, 2014 06:19 PM

A GeoSpatial World

Pg3DViewer nouvelle version

 

Pg3DViewer nouvelle version

image

J’ai réalisé une nouvelle version de l’outil Pg3DViewer qui intègre ScintillaNET un puissant contrôle d’édition de textes avec coloration syntaxique pour Windows et un wrapper pour le composant Scintilla. J’ai utilisé ce contrôle pour gérer la syntaxe SQL de PostgreSQL  et les fonctions de PostGIS. Le fichier pgsql.xml que vous trouverez dans le répertoire d’installation contient tous les éléments pour la gestion des mots clefs et du style affecté pour chaque liste de mots clefs, je vous laisse le plaisir de l’explorer et de le modifier selon vos désirs.

Le lien de téléchargement : Pg3DViewer nouvelle version

by Jérôme ROLLAND (noreply@blogger.com) at September 09, 2014 04:00 PM

September 02, 2014

Boston GIS (Regina Obe, Leo Hsu)

PostGIS sessions at Postgres Open 2014 Chicago not too late

If you haven't signed up for our Postgres Open 2014 PostGIS tutorials, there are still a few tickets available, but time and space is running out.

For those of you who've already signed up for our tutorials in September 17th, 2014, keep an eye on this page: http://www.postgis.us/pgopen2014. Not much to see there yet, but we'll be posting content there before the tutorials as well as the examples so you can flaunt your mastery of cutting and pasting.

by Regina Obe (nospam@example.com) at September 02, 2014 03:20 AM

Postgres OnLine Journal (Leo Hsu, Regina Obe)

PostGIS presentations at Postgres Open 2014 in Chicago

For those of you who will be attending either of our PostGIS sessions in Postgres Open 2014 in Chicago September 17th, we've put up a page where we'll be posting links to the data we'll be using as well as the examples so you can follow along and copy and paste. Page is http://www.postgis.us/pgopen2014.

We ask that you bookmark the page before you come and try to install PostGIS 2.1 if you haven't already. Not much to see there yet. For those who haven't signed up, it's not too late. For those who will be unable to make it, we'll try to post all the content from the tutorials on the above page.

by Leo Hsu and Regina Obe (nospam@example.com) at September 02, 2014 03:05 AM

August 20, 2014

shisaa.jp

Postgis and PostgreSQL in Action - Timezones

Preface

Recently, I was lucky to be part of an awesome project called the Breaking Boundaries Tour.

This project is about two brothers, Omar and Greg Colin, who take their Stella scooters to make a full round trip across the United States. And, while they are at it, try to raise funding for Surfer's Healing Folly Beach - an organization that does great work enhancing the lives of children with autism through surfing . To accommodate this trip, they wished to have a site where visitors could follow their trail live, as it happened. A marker would travel across the map, with them, 24/7.

Furthermore, they needed the ability to jump off their scooters, snap a few pictures, edit a video, write some side info and push it on the net, for whole the world to see. Immediately after they made their post, it had to appear on the exact spot they where at when snapping their moments of beauty.

To aid in the live tracking of their global position, they acquired a dedicated GPS tracking device which sends a latitude/longitude coordinate via a mobile data network every 5 minutes.

Now, this (short) post is not about how I build the entire application, but rather about how I used PostGIS and PostgreSQL for a rather peculiar matter: deducting timezone information.

For those who are interested though: the site is entirely build in Python using the Flask "micro framework" and, of course, PostgreSQL as the database.

Timezone information?

Yes. Time, dates, timezones: hairy worms in hairy cans which many developers hate to open, but have to sooner or later.

In the case of Breaking Boundaries Tour, we had one major occasion where we needed the correct timezone information: where did the post happen?

Where did it happen?

A feature we wanted to implement was one to help visitors get a better view of when a certain post was written. To be able to see when a post was written in your local timezone is much more convenient then seeing the post time in some foreign zone.

We are lazy and do not wish to count back- or forward to figure out when a post popped up in our frame of time.

The reasoning is simple, always calculate all the times involved back to simple UTC (GMT). Then figure out the clients timezone using JavaScript, apply the time difference and done!

Simple eh?

Correct, except for one small detail in the feature request, in what zone was the post actually made?

Well...damn.

While you heart might be at the right place while thinking: "Simple, just look at the locale of the machine (laptop, mobile phone, ...) that was used to post!", this information if just too fragile. Remember, the bothers are crossing the USA, riding through at least three major timezones. You can simply not expect all the devices involved when posting to always adjust their locale automatically depending on where they are.

We need a more robust solution. We need PostGIS.

But, how can a spatial database help us to figure out the timezone?

Well, thanks to the hard labor delivered to us by Eric Muller from efele.net, we have a complete and maintained shapefile of the entire world, containing polygons that represent the different timezones accompanied by the official timezone declarations.

This enables us to use the latitude and longitude information from the dedicated tracking device to pin point in which timezone they where while writing their post.

So let me take you on a short trip to show you how I used the above data in conjunction with PostGIS and PostgreSQL.

Getting the data

The first thing to do, obviously, is to download the shapefile data and load it in to our PostgreSQL database. Navigate to the Timezone World portion of the efele.net site and download the "tz_world" shapefile.

This will give you a zip which you can extract:

$ unzip tz_world.zip

Unzipping will create a directory called "world" in which you can find the needed shapefile package files.

Next you will need to make sure that your database is PostGIS ready. Connect to your desired database (let us call it bar) as a superuser:

$ psql -U postgres bar

And create the PostGIS extension:

CREATE EXTENSION postgis;

Now go back to your terminal and load the shapefile into your database using the original owner of the database (here called foo):

$ shp2pgsql -S -s 4326 -I tz_world | psql -U foo bar

As you might remember from the PostGIS series, this loads in the geometry from the shapefile using only simple geometry (not "MULTI..." types) with a SRID of 4326.

What have we got?

This will take a couple of seconds and will create one table and two indexes. If you describe your database (assuming you have not made any tables yourself):

public | geography_columns | view     | postgres
public | geometry_columns  | view     | postgres
public | raster_columns    | view     | postgres
public | raster_overviews  | view     | postgres
public | spatial_ref_sys   | table    | postgres
public | tz_world          | table    | foo
public | tz_world_gid_seq  | sequence | foo

You will see the standard PostGIS bookkeeping and you will find the tz_world table together with a gid sequence.

Let us describe the table:

\d tz_world

And get:

Column |          Type          |                       Modifiers                        
--------+------------------------+--------------------------------------------------------
gid    | integer                | not null default nextval('tz_world_gid_seq'::regclass)
tzid   | character varying(30)  | 
geom   | geometry(Polygon,4326) | 
Indexes:
    "tz_world_pkey" PRIMARY KEY, btree (gid)
    "tz_world_geom_gist" gist (geom)

So we have:

  • gid: an arbitrary id column
  • tzid: holding the standards compliant textual timezone identification
  • geom: holding polygons in SRID 4326.

Also notice we have two indexes made for us:

  • tz_world_pkey: a simple B-tree index on our gid
  • tz_world_geom_gist: a GiST index on our geometry

This is a rather nice set, would you not say?

Using the data

So how do we go about using this data?

As I have said above, we need to figure out in which polygon (timezone) a certain point resides.

Let us take an arbitrary point on the earth:

  • latitude: 35.362852
  • longitude: 140.196131

This is a spot in the Chiba prefecture, central Japan.

Using the Simple Features functions we have available in PostGIS, it is trivial to find out in which polygon a certain point resides:

SELECT tzid 
    FROM tz_world
    WHERE ST_Intersects(ST_GeomFromText('POINT(140.196131 35.362852)', 4326), geom);

And we get back:

    tzid    
------------
 Asia/Tokyo

Awesome!

In the above query I used the function ST_Intersects which checks if a given piece of geometry (our point) shares any space with another piece. If we would check the execute plan of this query:

EXPLAIN ANALYZE SELECT tzid 
    FROM tz_world
    WHERE ST_Intersects(ST_GeomFromText('POINT(140.196131 35.362852)', 4326), geom);

We get back:

                                                      QUERY PLAN                                                          
------------------------------------------------------------------------------------------------------------------------------
Index Scan using tz_world_geom_gist on tz_world  (cost=0.28..8.54 rows=1 width=15) (actual time=0.591..0.592 rows=1 loops=1)
    Index Cond: ('0101000020E61000006BD784B446866140E3A430EF71AE4140'::geometry && geom)
    Filter: _st_intersects('0101000020E61000006BD784B446866140E3A430EF71AE4140'::geometry, geom)
Total runtime: 0.617 ms

That is not bad at all, a runtime of little over 0.6 Milliseconds and it is using our GiST index.

But, if a lookup is using our GiST index, a small alarm bell should go off inside your head. Remember my last chapter on the PostGIS series? I kept on babbling about index usage and how geometry functions or operators can only use GiST indexes when they perform bounding box calculations.

The latter might pose a problem in our case, for bounding boxes are a very rough approximations of the actual geometry. This means that when we arrive near timezone borders, our calculations might just give us the wrong timezone.

So how can we fix this?

This time, we do not need to.

This is one of the few blessed functions that makes use of both an index and is very accurate.

The ST_Intersects first uses the index to perform bounding box calculations. This filters out the majority of available geometry. Then it performs a more expensive, but more accurate calculation (on a small subset) to check if the given point is really inside the returned matches.

We can thus simply use this function without any more magic...life is simple!

Implementation

Now it is fair to say that we do not wish to perform this calculation every time a user views a post, that would not be very efficient nor smart.

Rather, it is a good idea to generate this information at post time, and save it for later use.

The way I have setup to save this information is twofold:

  • I only save a UTC (GTM) generalized timestamp of when the post was made.
  • I made an extra column in my so-called "posts" table where I only save the string that represents the timezone (Asia/Tokyo in the above case).

This keeps the date/time information in the database naive of any timezone and makes for easier calculations to give the time in either the clients timezone or in the timezone the post was originally written. You simply have one "root" time which you can move around timezones.

On every insert of a new post I have created a trigger that fetches the timezone and inserts it into the designated column. You could also fetch the timezone and update the post record using Python, but opting for an in-database solution saves you a few extra, unneeded round trips and is most likely a lot faster.

Let us see how we could create such a trigger.

A trigger in PostgreSQL is an event you can set to fire when certain conditions are met. The event(s) that fire have to be encapsulated inside a PostgreSQL function. Let us thus first start by creating the function that will insert our timezone string.

Creating functions

In PostgreSQL you can write functions in either C, Procedural languages (PgSQL, Perl, Python) or plain SQL.

Creating functions with plain SQL is the most straightforward and most easy way. However, since we want to write a function that is to be used inside a trigger, we have even a better option. We could employ the power of the embedded PostgreSQL procedural language to easily access and manipulate our newly insert data.

First, let us see which query we would use to fetch the timezone and update our post record:

UPDATE posts
  SET tzid = timezone.tzid
  FROM (SELECT tzid
          FROM tz_world
          WHERE ST_Intersects(
            ST_SetSRID(
              ST_MakePoint(140.196131, 35.362852),
            4326),
          geom)) AS timezone
  WHERE pid = 1;

This query will fetch the timezone string using a subquery and then update the correct record (a post with "pid" 1 in this example).

How do we pour this into a function?

CREATE OR REPLACE FUNCTION set_timezone() RETURNS TRIGGER AS $$
BEGIN
    UPDATE posts
    SET tzid = timezone.tzid
    FROM (SELECT tzid 
            FROM tz_world
              WHERE ST_Intersects(
                ST_SetSRID(
                  ST_MakePoint(NEW.longitude, NEW.latitude),
                4326),
              geom)) AS timezone
    WHERE pid = NEW.pid;
    RETURN NEW;
END $$
LANGUAGE PLPGSQL IMMUTABLE;

First we use the syntax CREATE OR REPLACE FUNCTION to indicate we want to create (or replace) a custom function. Then we tell PostgreSQL that this function will return type TRIGGER.

You might notice that we do not give this function any arguments. The reasoning here is that this function is "special". Functions which are used as triggers magically get information about the inserted data available.

Inside the function you can see we access our latitude and longitude prefixed with NEW. These keywords, NEW and OLD, refer to the record after and before the trigger(s) happened. In our case we could have used both, since we do not alter the latitude or longitude data, we simply fill a column that is NULL by default. There are more keywords available (TG_NAME, TG_RELID, TG_NARGS, ...) which refer to properties of the trigger itself, but that is beyond today's scope.

The actual SQL statement is wrapped between double dollar signs ($$). This is called dollar quoting and is the preferred way to quote your SQL string (as opposed to using single quotes). The body of the function, which in our case is mostly the SQL statement, is surrounded with a BEGIN and END keyword.

A trigger function always needs a RETURN statement that is used to provide the data for the updated record. This too has to reside in the body of the function.

Near the end of our function we need to declare in which language this function was written, in our case PLPGSQL.

Finally, the IMMUTABLE keyword tells PostgreSQL that this function is rather "functional", meaning: if the inputs are the same, the output will also, always be the same. Using this caching keyword gives our famous PostgreSQL planner the ability to make decisions based on this knowledge.

Creating triggers

Now that we have this functionality wrapped into a tiny PLPGSQL function, we can go ahead and create the trigger.

First you have the event on which a trigger can execute, these are:

  • INSERT
  • UPDATE
  • DELETE
  • TRUNCATE

Next, for each event you can specify at what timing your trigger has to fire:

  • BEFORE
  • AFTER
  • INSTEAD OF

The last one is a special timing by which you can replace the default behavior of the mentioned events.

For our use case, we are interested in executing our function AFTER INSERT.

CREATE TRIGGER set_timezone
    AFTER INSERT ON posts
    FOR EACH ROW
    EXECUTE PROCEDURE set_timezone();

This will setup the trigger that fires after the insert of a new record.

Wrapping it up

Good, that all there is to it.

We use a query, wrapped in a function, triggered by an insert event to inject the official timezone string which is deducted by PostGIS's spatial abilities.

Now you can use this information to get the exact timezone of where the post was made and use this to present the surfing client both the post timezone time and their local time.

For the curious ones out there: I used the MomentJS library for the client side time parsing. This library offers a timezone extension which accepts these official timezone strings to calculate offsets. A lifesaver, so go check it out.

Also, be sure to follow the bros while they scooter across the States!

And as always...thanks for reading!

by Tim van der Linden at August 20, 2014 10:00 AM

August 08, 2014

Stephen Mather

KNN with FLANN and laspy, a starting place

FLANN is Fast Library for Approximate Nearest Neighbors, which is a purportedly wicked fast nearest neighbor library for comparing multi-dimensional points. I only say purportedly, as I haven’t verified, but I assume this to be quite true. I’d like to move some (all) of my KNN calculations outside the database.

I’d like to do the following with FLANN– take a LiDAR point cloud and change it into a LiDAR height-above-ground point cloud. What follows is my explorations so far.

In a previous series of posts, e.g. http://smathermather.wordpress.com/2014/07/14/lidar-and-pointcloud-extension-pt-6/

pointcloud6

I have been using the point cloud extension in PostGIS. I like the 2-D chipping, but I think I should segregate my data into height classes before sending it into the database. In this way, I can query my data by height class and by location efficiently, taking full advantage of the efficiencies of storing all those little points in chips, while also being able to query the data in any of the dimensions I need to in the future. Enter FLANN.

I haven’t gotten far. To use FLANN with LiDAR through Python, I’m also using laspy.  There’s a great tutorial here: http://laspy.readthedocs.org/en/latest/tut_part_1.html

laspy

I make one change to the tutorial section using FLANN. The code as written is:

import laspy
import pyflann as pf
import numpy as np

# Open a file in read mode:
inFile = laspy.file.File("./laspytest/data/simple.las")
# Grab a numpy dataset of our clustering dimensions:
dataset = np.vstack([inFile.X, inFile.Y, inFile.Z]).transpose()

# Find the nearest 5 neighbors of point 100.

neighbors = flann.nn(dataset, dataset[100,], num_neighbors = 5)
print("Five nearest neighbors of point 100: ")
print(neighbors[0])
print("Distances: ")
print(neighbors[1])

To make this example work with the current version of pyflann, we need to make sure we import all of pyflann (or at least nn), and also set flann = FLANN() as follows:

import laspy
import numpy as np
from pyflann import *

# Open a file in read mode:
inFile = laspy.file.File("simple.las")
# Grab a numpy dataset of our clustering dimensions:
dataset = np.vstack([inFile.X, inFile.Y, inFile.Z]).transpose()

# Find the nearest 5 neighbors of point 100.
flann = FLANN()
neighbors = flann.nn(dataset, dataset[100,], num_neighbors = 5)
print("Five nearest neighbors of point 100: ")
print(neighbors[0])
print("Distances: ")
print(neighbors[1])

Finally, a small note on installation of pyflann on Ubuntu. What I’m about to document is undoubtedly not the recommended way to get pyflann working. But it worked… .

Installation for FLANN on Ubuntu can be found here: http://www.pointclouds.org/downloads/linux.html
pcl

But this does not seem to install pyflann. That said, it installs all our dependencies + FLANN, so…

I cloned, compiled, and installed the FLANN repo: https://github.com/mariusmuja/flann

git clone git://github.com/mariusmuja/flann.git
cd flann
mkdir BUILD
cd BUILD
cmake ../.
make
sudo make install

This get’s pyflann where it needs to go, and voila! we can now do nearest neighbor searches within Python.

Next step, turn my LiDAR xyz point cloud into a xy-height point cloud, then dump in height-class by height class into PostgreSQL. Wish me luck!


by smathermather at August 08, 2014 02:15 PM

August 07, 2014

Stephen Mather

Drivetime analyses, pgRouting

Map of overlapping 5-minute drive times from park access points
We’ve got some quick and dirty pgRouting-based code up on github. I say quick and dirty because it directly references the table names in both of the functions. I hope to fix this in the future.

The objective with this code is to input a point, use a nearest neighbor search to find the nearest intersection, and from that calculate a drive time alpha shape.

First, the nearest neighbor search:

CREATE OR REPLACE FUNCTION pgr.nn_search (geom geometry) RETURNS
int8 AS $$

SELECT id FROM pgr.roads_500_vertices_pgr AS r
ORDER BY geom <#> r.the_geom
LIMIT 1;

$$ LANGUAGE SQL VOLATILE;

This function takes one argument– pass it a point geometry, and it will do a knn search for the nearest point. Since we are leveraging pgRouting, it simply returns the id of the nearest intersection point. We wrote this as a function in order to run it, e.g. on all the points in a table. As I stated earlier, it directly references a table name, which is a little hackerish, but we’ll patch that faux pax later.

Now that we have the ID of the point in question, we can do a driving distance calculation, wrap that up in an alpha shape, and return our polygon. Again, we write this as a function:

CREATE OR REPLACE FUNCTION pgr.alpha_shape (id integer, minutes integer) RETURNS
geometry AS $$

WITH alphashape AS(SELECT pgr_alphaShape('WITH
DD AS (
SELECT seq, id1 AS node, cost
FROM pgr_drivingDistance(''SELECT id, source, target, cost FROM pgr.roads_500'',' || id || ', ' || minutes || ', false, false)),
dd_points AS (
SELECT id_ AS id, x, y
FROM pgr.roads_500_vertices_pgr v, DD d
WHERE v.id = d.node)
SELECT * FROM dd_points')),

alphapoints AS (
SELECT ST_Makepoint((pgr_alphashape).x, (pgr_alphashape).y) FROM alphashape),

alphaline AS (
SELECT ST_MakeLine(ST_MakePoint) FROM alphapoints)

SELECT ST_MakePolygon(ST_AddPoint(ST_MakeLine, ST_StartPoint(ST_MakeLine))) AS the_geom FROM alphaline

$$ LANGUAGE SQL VOLATILE;

Finally, we’ll use these functions in conjunction with a set of park feature access points to map our our 5-minute drive time. Overlaps in 5-minute zones we’ve rendered as a brighter green in the image above.

CREATE TABLE pgr.alpha_test AS
WITH dest_ids AS (
SELECT pgr.nn_search(the_geom) AS id FROM pgr.dest_pts
)
SELECT pgr.alpha_shape(id::int, 5)::geometry, id FROM dest_ids;

Oh, and I should point out that for the driving distance calculation, we have pre-calculated the costs of the roads based on the speed limit, so cost is really travel time in minutes.


by smathermather at August 07, 2014 01:50 AM

August 01, 2014

Stephen Mather

Plugin-free QGIS TMS tiles via GDAL

Want to load your favorite tiles into QGIS? How about a plugin-free QGIS TMS tiles via GDAL:

http://www.3liz.com/blog/rldhont/index.php?post/2012/07/17/OpenStreetMap-Tiles-in-QGIS

Really awesome… .

Needs but one change: epsg:900913 should be epsg:3857 or QGIS (GDAL?) throws an error. Presumably you could also define epsg:900913 in some config file, but barring that create an XML file as follows, and load as a raster in QGIS:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://maps1.clemetparks.com/tilestache/tilestache.cgi/basemap/${z}/${x}/${y}.jpg</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>

Now I can use my pretty tilestache maps _everywhere!_:
Image of coyote points overlayed on tilestache rendered hillshade map

Screenshot of map from postgis

edit: forgot the hattip for the lead on this, Thomas gratier: @ThomasG77


by smathermather at August 01, 2014 03:28 AM

Stephen Mather

Cleaning animal tracking data — throwing away extra points

Much the problem of the modern era– too much data, uneven data, and yet, should we keep it all?

Here’s the problem space: attach GPS collar to a coyote, send that data home, and you have a recipe for gleaning a lot of information about the movement of that animal across the landscape. In order to maximize the data collected while also maximizing the battery life of said device, sometimes (according to a preset scheme) you might sample every hour or every 15 minutes, and then switch back to once a day.

When you get the data back in, you get something like this:

Map of raw coyote data-- all points, showing all available collar GPS points

Already, we see some pretty interesting patterns. The first thing that I notice is clusters of activity. Are these clusters related to our uneven sampling, sometimes every 15 minutes, sometimes once or twice a day? Or are those really important clusters in the overall activity of the animal?

The next thing I notice is how perfectly the range of movement of the coyote is bounded by expressways to the west and north, and by the river to the east. Those seem to be pretty impermiable boundaries.

So, as does a man who only has a hammer, we clean the data with PostGIS. In this case, it’s an ideal task. First, metacode:

Task one: union / dissolve our points by day — for some days this will give us a single point, for other days a multi-point for the day.

Task two: find the centroid of our union. This will grab a centroid for each clump of points (a weighted spatial mean, if you will) and return that single point per day. There are some assumptions here that might not bear out under scrutiny, but it’s not a bad first approximation.

Now, code:

CREATE TABLE centroids AS
-- union / dissolve our points by day
WITH clusterp AS (
	SELECT ST_Union(geom) AS geom, gmt_date FROM nc_clipped_data_3734
		GROUP BY gmt_date
),
-- find the centroid of our union
centroids AS (
	SELECT ST_Centroid(geom) AS geom, gmt_date FROM clusterp
)
SELECT * FROM centroids;

Boom! One record per day:

Map showing much smaller home range, one record per day clustered around a single core zone

A different pattern now emerges. We can’t see the hard boundaries of the use area, but now the core portion of the home range can be seen.

Upcoming (guest) blog post on the use of these data in home range calculations in R. Stay tuned.

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA.


by smathermather at August 01, 2014 02:58 AM

July 28, 2014

A GeoSpatial World

PostGIS Viewer 3D suite

 

PosGIS Viewer 3D suite

  pg3DViewer_title6_thumb8

image_thumb175 ntroduction

Comme je vous l’annoncez dans mon précédent billet j’ai fait évoluer le pg3DViewer pour qu’il puisse gérer les géométries de type POINT, et pour mes test je me suis intéressé aux nuages de points ‘Points Clouds’'

image

 

 

 

 

 

 

 

Pour la suite de ce tutorial nous utiliserons deux fichiers de données dont voici les liens de téléchargement :

image_thumb174 nstallation

image_thumb41

Cliquez sur le lien pg3DViewer_setup.exe pour télécharger le programme d’installation de la nouvelle version, si vous aviez installé la version précédente désinstallez la, puis lancer l’exécution.

Cliquer sur le  bouton Suivant à chaque étape, puis sur le bouton Terminer.

Dans le précédent billet vous avez tous les écrans de l’installation : http://ageoguy.blogspot.fr/2014/07/postgis-3d-viewer.html

 

image_thumb179 tilisation

Avant de commencer à utiliser pg3Dviewer vous devez avoir une base de données PostgreSQL avec la cartouche spatiale PostGIS version 2.0 ou supérieure.

image_thumb61 Double cliquez sur l’icone créé par le programme d’installation sur le bureau pour lancer l’application.

image_thumb50

image_thumb2

Connecter vous à votre base de données 3D

 

image_thumb16

Commencez par cliquer sur cet icône pour vous connecter a une base de données PostgreSQL/PostGIS version 2.0 ou supérieur contenant des données 3D ou prête à en recevoir.

 
  • Choisissez localhost pour serveur si votre base de données est en local sur votre machine ou bien saisissez l’adresse IP de votre serveur
  • le port est par défaut 5432, changez le si nécessaire
  • saisissez l’utilisateur
  • saisissez le mot de passe
  • choisissez une base de donnée dans la combobox
  • cliquez sur le bout OK pour valider votre saisie.

image_thumb1711

 

Traitement d’un fichier Ascii XYZ

Nous allons utiliser le fichier bunny.xyz que vous avez téléchargé dans l’introduction.

Nous allons créer la table qui va accueillir les donnés. Dans PgAdmin III ouvrez l’outil de requêtes puis
copiez la requête suivante :

-- Table: bunny

-- DROP TABLE bunny;

CREATE TABLE bunny
(
  x double precision,
  y double precision,
  z double precision,
  id bigserial NOT NULL,
  CONSTRAINT bunny_pkey PRIMARY KEY (id)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE bunny
  OWNER TO postgres;

Puis exécuter la.

image
Puis nous allons insérer les données dans la table. Copiez la requête suivante dans l’éditeur de requêtes :

copy lidar_test(x,y,z,classification)
from 'C:/temp/bunny.txt'
with delimiter ' '

puis exécuter la.

image
Maintenant nous allons pouvoir visualiser les données. Dans le panneau Query copiez la requête suivante :

SELECT x,y,z
FROM bunny AS points_xyz

Cliquez sur l’icone ‘Execute query

image
Et voici le résultat après Zoom et Rotation.
Comme vous avez du le remarquer le temps d’affichage est très rapide pour une table contenant 34835 points.














image
Nous allons maintenant essayer une nouvelle requête.Dans le panneau Query copiez la requête suivante :

SELECT st_geomfromtext('POINT Z('||x||' '||y||' '||z||')',0) as geom,'255:0:0'::text AS rgb FROM bunny

Cliquez sur l’icone ‘Execute query

image
Et voici le résultat après Zoom et Rotation.
Le temps d’affichage est beaucoup plus long parce que dans cette requête nous construisons la géométrie à la volée.

La différence avec la requête précédente est que celle-ci génère un fichier de points au format xyz (grâce à la commande COPY TO de PostgreSQL) qu’exploite un ‘reader’ de  la bibliothèque VTK. Tout cela ce fait avec l’ajout de ‘A'S points_xyz’ a la fin de la requête et la sélection des champs x,y et z.





image

Traitement d’un fichier LAS

Nous allons utiliser le fichier ‘2013 LiDAR 5105_54480.las’ que vous avez téléchargé dans l’introduction. Il va nous falloir des outils pour pouvoir importer ce fichier dans notre base de donnée, nous allons utiliser une collection d’outils qui s’appelle LASTOOLS que l’on peut télécharger sur le site de http://www.cs.unc.edu/~isenburg/lastools/ et dont voici le lien de téléchargement : http://lastools.org/download/lastools.zip

L’outil que nous allons utiliser s’appelle las2txt , voici l’entête du README associé:

****************************************************************
las2txt:
Converts from binary LAS/LAZ 1.0 - 1.4 to an ASCII text format
For updates check the website or join the LAStools mailing list.

http://rapidlasso.com/
http://lastools.org/
http://groups.google.com/group/lastools/
http://twitter.com/lastools/
http://facebook.com/lastools/
http://linkedin.com/groups?gid=4408378

Martin @lastools
****************************************************************

Et le lien pour le consulter : http://www.cs.unc.edu/~isenburg/lastools/download/las2txt_README.txt

Ouvrez une fenêtre DOS, allez dans le répertoire ou vous avez installé LASTOOLS, puis dans le sous répertoire bin. Copiez la commande suivante :
las2txt -i "C:\temp\2013 LiDAR 5105_54480.las" -parse xyzc -sep semicolon

Cette commande va générer un fichier Ascii (.txt) au format x,y,z et classification avec pour séparateur de colonne le caractère ‘;’ (point virgule).

image
Nous disposons à présent d’un fichier que nous allons pouvoir insérer dans notre base de données. Pour cela nous allons commencer par créer la table qui recevra les données, dans la fenêtre de requêtes de PgAdmin III copiez la requête suivante:

-- Table: lidar_test

-- DROP TABLE lidar_test;

CREATE TABLE lidar_test
(
  x double precision,
  y double precision,
  z double precision,
  classification integer
)
WITH (
  OIDS=FALSE
);
ALTER TABLE lidar_test
  OWNER TO postgres;

puis exécuter la.
image

Puis nous allons insérer les données dans la table. Copiez la requête suivante dans l’éditeur de requêtes :

copy lidar_test(x,y,z,classification)
from 'C:/temp/2013 LiDAR 5105_54480.txt'
with delimiter ';'

puis exécuter la.

image
Maintenant nous allons pouvoir visualiser les données. Dans le panneau Query copiez la requête suivante :

select x,y,z,
CASE
WHEN classification=2 then '128:0:0'
WHEN classification=3 then '0:64:0'
WHEN classification=4 then '0:128:0'
WHEN classification=5 then '0:255:0'
WHEN classification=6 then '255:128:64'
WHEN classification=11 then '128:128:128'
END as rgb
from lidar_test as points_xyzrgb

Cliquez sur l’icone ‘Execute query

image
Et voici le résultat
















image
Les valeurs de classification utilisées dans la requête sont issues de cette requête :
select distinct classification from lidar_test order by 1

Vous pouvez également obtenir des informations sur le fichier LAS avec l’outil lasinfo, dans une fenêtre DOS tapez la commande suivante :

lasinfo c:\temp\2013 LiDAR 5105_54480.las

a la fin du rapport vous obtenez la définition de la classification.

lasinfo for c:\temp\2013 LiDAR 5105_54480.las
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            1
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          ''
  generating software:        'TerraScan'
  file creation day/year:     156/2013
  header size:                227
  offset to point data:       447
  number var. length records: 2
  point data format:          1
  point data record length:   28
  number of point records:    7938187
  number of points by return: 6742290 914660 231891 44542 4804
  scale factor x y z:         0.001 0.001 0.001
  offset x y z:               510499.94 5447999.8799999999 0
  min x y z:                  510499.940 5447999.880 34.110
  max x y z:                  511000.110 5448500.100 275.690
variable length header record 1 of 2:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  88
  description          'Georeferencing Information'
    GeoKeyDirectoryTag version 1.1.0 number of keys 10
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 2048 tiff_tag_location 0 count 1 value_offset 4269 - GeographicTypeGeoKey: GCS_NAD83
      key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
      key 2056 tiff_tag_location 0 count 1 value_offset 7019 - GeogEllipsoidGeoKey: Ellipse_GRS_1980
      key 2057 tiff_tag_location 34736 count 1 value_offset 0 - GeogSemiMajorAxisGeoKey: 6378137
      key 2058 tiff_tag_location 34736 count 1 value_offset 1 - GeogSemiMinorAxisGeoKey: 6356752.314
      key 2059 tiff_tag_location 34736 count 1 value_offset 2 - GeogInvFlatteningGeoKey: 298.2572221
      key 3072 tiff_tag_location 0 count 1 value_offset 26910 - ProjectedCSTypeGeoKey: UTM 10 northern hemisphere
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
variable length header record 2 of 2:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34736
  length after header  24
  description          'Double Param Array'
    GeoDoubleParamsTag (number of doubles 3)
      6.37814e+006 6.35675e+006 298.257
reporting minimum and maximum for all LAS point record entries ...
  X                   0     500170
  Y                   0     500220
  Z               34110     275690
  intensity           0        255
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          0
  scan_direction_flag 0          1
  classification      2         11
  scan_angle_rank   -19         22
  user_data          79         79
  point_source_ID 10217      10238
  gps_time 49786940.848411 49796462.877692
WARNING: there is coordinate resolution fluff (x10) in XYZ
overview over number of returns of given pulse: 5825940 1363444 563503 160908 24392 0 0
histogram of classification of points:
         2030423  ground (2)
         1836807  low vegetation (3)
          234836  medium vegetation (4)
         1811842  high vegetation (5)
         2022021  building (6)
            2258  road surface (11)

 

image_thumb100 onclusion.

Nous sommes arrivés à la fin de ce billet, je vous conseille la lecture d’un excellent tutorial de Paul Ramsey concernant les données LIDAR que vous trouverez a cette adresse  http://workshops.boundlessgeo.com/tutorial-lidar/ , il ne vous restera plus qu’a afficher les données en utilisant les fonctions présentées dans ce ‘tutorial lidar’.

by Jérôme ROLLAND (noreply@blogger.com) at July 28, 2014 10:08 AM

July 25, 2014

A GeoSpatial World

Plugin PgAdmin III : PostGISViewer suite

Plugin PgAdmin III : PostGISViewer  suite

Nouvelle version

 PostGISViewer

Multi-géométrie

J'ai mis en place la prise en compte des multi-géométries  :
  • MULTIPOLYGON
  • MULTILINESTRING
  • MULTIPOINT

Requêtes
  
J'ai rajouté la possibilité d'exécuter les requêtes présente dans l'éditeur SQL de PgAdmin III, mais pour que cela soit possible, il faut respecter certaines règles:
  • Un seul éditeur SQL de PgAdmin III doit être ouvert.
  • Pour que la ou les requête(s) soient exécutées, il ne faut pas être positionné sur un nom de table dans le navigateur d'objets de PgAdmin III , ou alors que la table soit déjà chargée dans le visualiseur.  Une table sélectionnée prendra toujours le pas sur une ou des requêtes.
  • Le premier champ doit être  géométrique ou une géométrie issue d'une fonction de traitement géométrique (ST_Buffer, ST_Intersection....) par requête sera utilisé pour l'affichage. Ce champ ou cette géométrie sera encapsulé par la fonction AsBinary et aura pour alias geom :
    • SELECT AsBinary(wkb_geometry) as geom ....
    • SELECT AsBinary(ST_Buffer(wkb_geometry, 1.0)) as geom...
    • SELECT AsBinary(ST_Intersection(a.wkb_geometry,b.wkb_geoemetry)) as geom... 
    • ....
  •  Le second champ doit permettre de déterminer le type de géométrie, il doit donc être encapsulé par la fonction GeometryType, un alias n'est pas nécessaire :
    • SELECT ..., GeometryType(wkb_geometry)...
    • SELECT ..., GeometryType(ST_Buffer(wkb_geometry,1.0))...
    • SELECT ..., GeometryType(ST_Intersection(a.wkb_geometry,b.wkb_geometry))...
    • ....

  • Chaque requête devra se terminer par un point virgule, ce qui permettra de pouvoir exécuter plusieurs requêtes a la suite.
 Exemple

Ci-dessous deux requêtes se terminant par des points virgules :
  • La première requête va charger la commune qui a pour nom 'Sainte-Croix-de-Quintillargues'
  • La seconde requête va charger tous les bâtiments de cette commune.
 

































Après avoir lancé le plugin PostGISViewer, les deux requêtes sont exécutées et donne le résultat suivant  :





































Les couches créées portent comme nom query avec un identifiant qui s'incrémente pour toute la cession du visualiseur.

Toutes suggestions, remarques pour l'amélioration de cet outil seront les bienvenues.

A suivre...

by Jérôme ROLLAND (noreply@blogger.com) at July 25, 2014 02:48 PM

July 20, 2014

A GeoSpatial World

PostGIS 3D Viewer

 

image n viewer 3D pour PostGIS

image pg3DViewer_title[6]
 

image ntroduction

Avec la version 2.0 PostGIS est entré dans le monde de la 3D, il est possible de stocker des géométries 3D avec l’ajout de la coordonnée Z dans toutes les géométries existantes et de stocker des géométries volumétriques avec l’introduction de deux nouveaux type  :

  • TIN(Triangulated Irregular Network, une collection de TRIANGLE)
image
  • POLYHEDRALSURFACE (une collection de POLYGON)
image

Cela fait déjà un an que la création d’un viewer pour visualiser les géométries 3D de PostGIS me trottait dans la tête suite à un échange avec Claus Nagel et Felix Kunde de virtualcitySYSTEMS à propos de l’outil d’import du format CityGML que j’avais réalisé CityGML vers PostGIS 2.0 , ils m’avaient conseillé alors de me diriger vers le développement d’un viewer 3D. Lorsque récemment j’ai repris contact avec eux pour leur parler de mon développement en cours, ils ont accueilli la nouvelle avec enthousiasme et ont bénéficié de versions beta, ce qui m’a grandement aidé pour mes choix dans la finalisation de mon application.

Aujourd’hui la 3D est de plus en plus présente dans le monde SIG avec le format CityGML et dans le monde de l’architecture avec le format IFC, importer ces formats dans une base de données PosgtgreSQL/PostGIS est une démarche naturelle car cela apporte toute la puissance de traitement des bases de données spatiale pour la gestion de ces données.

Pour le format CityGML la société virtualcitySYSTEMS met à disposition un schéma de base de données Open Source et gratuit pour le stockage et la gestion de ‘3D City Models’. Ce schéma de base de données est le reflet du modèle de données CityGML. Les données peuvent être importées via un outil d’import/export .
le lien pour accéder à leur site : http://virtualcitysystems.de/en/solutions.html#3dcitydb 
image

A ce stade du développement je n’affiche que des géométries de type POLYGON, MULTIPOLYGON, TIN et POLYHEDRALSURFACE. Ce qui permet d’afficher toutes type de géométries volumétriques répondant à cette contrainte. Je vais travailler rapidement sur les autres type de géométrie POINT, LINESTRING… , ce qui fera l’objet d’une nouvelle version de l’application dans les prochaines semaines.

J’ai développé cet outil en C# 2013 avec :

  • La librairie GDAL/OGR
[image66.png] GDAL - Geospatial Data Abstraction Library
  • ActiViz .Net un wrapper pour VTK
  ActiViz .Net


vtk-logo VTK
  • PostgreSQL 9.2
(thumbnail) PostgreSQL Windows
  • PostGIS 2.2.0 dev
adbadge_wide_240 Postgis
http://winnie.postgis.net/download/windows/pg92/buildbot/postgis-pg92-binaries-2.2.0devw32.zip
http://winnie.postgis.net/download/windows/pg92/buildbot/postgis-pg92-binaries-2.2.0devw64.zip

 

Pour l’intégration de données 3D j’ai testé l’outil 3DCityDB-Importer-Exporter http://http://www.3dcitydb.org/3dcitydb/downloads/ qui permet d’importer des fichiers au format CityGML dans une base de données PostgreSQL avec PostGIS (a partir de la version 2.0).

site_3dcitydb_02 

Voici une extraction du document 3DCityDB-v2_0_6-postgis-Tutorial.pdf concernant la partie PostGIS (en anglais ) :

image

 

image

image

image

image

image

image

image

image

 

image nstallation

image

Cliquez sur le lien pg3DViewer_setup.exe pour télécharger le programme d’installation, puis lancer l’exécution.

Cliquer sur le  bouton Suivant à chaque étape, puis sur le bouton Terminer.

pg3DViewer_install_01

pg3DViewer_install_02

pg3DViewer_install_03

pg3DViewer_install_04

pg3DViewer_install_05

pg3DViewer_install_06

 

image tilisation

Avant de commencer à utiliser pg3Dviewer vous devez avoir une base de données PostgreSQL avec la cartouche spatiale PostGIS version 2.0 ou supérieure. Cette base de données doit contenir des données 3D, l’idéal serait que vous ayez créé une base 3DCityDB et importé un ou plusieurs fichiers CityGML. Voici un lien ou vous pourrez télécharger des données CityGML Rotterdam 3D, demandez à notre ami Google de traduire la page cela peut aider, cela permet de découvrir un lien sur le site permettant de télécharger un fichier pdf contenant une carte avec tous les nom des quartiers téléchargeables http://www.rotterdam.nl/GW/Images/Gemeentewerken%20Home/Landmeten/overzicht%20buurten.pdf , il ne vous reste plus qu’a choisir les quartiers dans la liste sous le lien.

Vous pouvez quand même faire les premiers tests sans avoir téléchargé de données.

image Double cliquez sur l’icone créé par le programme d’installation sur le bureau pour lancer l’application.

image

image

Connecter vous à votre base de données 3D

 

image

Commencez par cliquer sur cet icône pour vous connecter a une base de données PostgreSQL/PostGIS version 2.0 ou supérieur contenant des données 3D ou prête à en recevoir.

 
  • Choisissez localhost pour serveur si votre base de données est en local sur votre machine ou bien saisissez l’adresse IP de votre serveur
  • le port est par défaut 5432, changez le si nécessaire
  • saisissez l’utilisateur
  • saisissez le mot de passe
  • choisissez une base de donnée dans la combobox
  • cliquez sur le bout OK pour valider votre saisie.

image


Visualiser une géométrie de type POLYHEDRALSURFACE


Dans le panneau Query copiez la requête suivante :

SELECT st_force_collection(
    st_geomFromText
    (
      
'POLYHEDRALSURFACE(
        ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
        ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
        ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
        ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)),
        ((1 0 0, 1 1 0, 1 1 1, 1 0 1, 1 0 0)),
        ((1 1 0, 0 1 0, 0 1 1, 1 1 1, 1 1 0))
        )'
    )
)

image

Cliquez sur l’icone ‘Execute query

image
Vous obtenez ce résultat.

Interactions de la souris avec la fenêtre de rendu:
- Bouton gauche : Rotation
- Roulette          : Zoom
- Bouton droit    : Zoom dynamique
- Bouton gauche + touche majuscule : Pan
Testez les différentes associations entre les touches Majuscule, Ctrl et alt et  les boutons de la souris.

Interactions touches claviers avec la fenêtre de rendu :
- Touche w : passe en mode affichage fil de fer.
- Touche s  : passe en mode affichage surfacique.
image
La fenêtre de rendu après rotation. image


Avant de continuer une petite explication sur la requête utilisé précédemment :

SELECT st_force_collection(
    st_geomFromText
    (
       
'POLYHEDRALSURFACE(
        ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
        ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
        ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
        ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)),
        ((1 0 0, 1 1 0, 1 1 1, 1 0 1, 1 0 0)),
        ((1 1 0, 0 1 0, 0 1 1, 1 1 1, 1 1 0))
        )'

    )
)

Comme vous avez du le remarquer la fonction ST_GeomFromText est encadré par la fonction ST_Force_Collection. Pourquoi rajouter cette fonction. Si vous exécuter cette requête sous pgAdmin III en encadrant le tout par ST_AsText vous obtenez comme résultat ceci :
"GEOMETRYCOLLECTION Z (
POLYGON Z ((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),
POLYGON Z ((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),
POLYGON Z ((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),
POLYGON Z ((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),
POLYGON Z ((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),
POLYGON Z ((1 1 0,0 1 0,0 1 1,1 1 1,1 1 0)))
"
Nous sommes passé d’une géométrie de type POLYHEDRALSURFACE à une géométrie de type GEOMETRYCOLLECTION Z contenant des POLYGON Z, mais rappelez vous une géométrie de type POLYHEDRALSURFACE est aussi une collection de POLYGON donc jusque la tout va bien.
Revenons à la question initiale pourquoi utiliser la fonction ST_Force_Collection, tout simplement parce la librairie GDAL que j’utilise dans mon développement ne reconnait pas encore les géométries de type POLYHDERALSURFACE. Cette astuce permet de contourner la limitation présente.

Visualiser une géométrie de type TIN

Dans le panneau Query copiez la requête suivante :

WITH res AS
(
SELECT st_polygon(
    st_boundary(
        (st_dump(
            st_geomFromText
            (
          
'TIN(
            ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
            ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
            ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
            ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'

            )
        )).geom
    )
    ,0
    ) as geom
)
SELECT st_forceCollection(st_collect(geom))
FROM res

image

Cliquez sur l’icone ‘Execute query

Vous obtenez ce résultat après rotation.

image


La requête utilisée nécessite une explication :

WITH res AS
(
SELECT st_polygon(
    st_boundary(
        (st_dump(
            st_geomFromText
            (
          
'TIN(
            ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
            ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
            ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
            ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'

            )
        )).geom
    )
    ,0
    ) as geom
)
SELECT st_forceCollection(st_collect(geom))
FROM res

Comme pour la requête concernant une géométrie de type POLYHEDRALSURFACE, les géométries de type TIN ne sont pas encore reconnues par la librairie GDAL. Pour rappel une géométrie de type TIN est aussi une collection de TRIANGLE que ne reconnait pas la librairie GDAL ce qui explique la requête légèrement plus complexe.
Si vous exécutez la requête suivante sous PgAdmin III :

SELECT St_AsText(st_force_collection(
    st_geomFromText
    (
         'TIN(
            ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
            ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
            ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
            ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )
'
    ))
)

Vous obtenez ceci :

"GEOMETRYCOLLECTION Z (
TRIANGLE Z ((0 0 0,1 0 0,0 1 0,0 0 0)),
TRIANGLE Z ((0 0 0,1 0 0,0 0 1,0 0 0)),
TRIANGLE Z ((0 0 0,0 0 1,0 1 0,0 0 0)),
TRIANGLE Z ((0 0 1,1 0 0,0 1 0,0 0 1)))
"

Que ne reconnait pas la librairie GDAL, ce qui explique la première partie de la requête qui recrée des géométries de type POLYGON à partir des géométries de type TRIANGLE. La seconde partie de la requête (ou l’on retrouve la fonction ST_Force_Collection) recrée une GEOMETRYCOLLECTION Z contenant des POLYGON Z donc interprétable par la librairie GDAL.

Pour simplifier les requêtes sur des géométries de type TIN, la création d’une fonction opérant la transformation me parait indispensable.

 

Fonction passant d’une géométrie de type TIN à une géométrie de type GEOMETRYCOLLECTION Z

-- Function: ST_TinToGeomCollection(geometry, INTEGER)

-- DROP FUNCTION ST_TinToGeomCollection(geometry, INTEGER);

CREATE OR REPLACE FUNCTION ST_TinToGeomCollection(geometry, integer)
  RETURNS geometry AS
$BODY$
DECLARE
    sql text;
    rec record;
BEGIN 
    sql:=
'WITH res AS
         (
          SELECT st_polygon(st_boundary((st_dump(geometry('||quote_literal($1::text)||'))).geom),'||$2||') as geom
          )
         SELECT st_forceCollection(st_collect(geom)) as geom FROM res
';
    FOR rec IN EXECUTE sql LOOP
       RETURN rec.geom;
   
END LOOP;
END
;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100;
ALTER FUNCTION ST_TinToGeomCollection(geometry,integer)
  OWNER TO postgres;
GRANT EXECUTE ON FUNCTION ST_TinToGeomCollection(geometry, integer) TO public;
GRANT EXECUTE ON FUNCTION ST_TinToGeomCollection(geometry, integer) TO postgres;

Vous allez créer la fonction dans votre base de données 3D :

  • Copier le contenu du cadre ci-dessus et collez dans la fenêtre de requêtes de PgAdmin III
  • lancez l’exécution, la fonction est créée.

la requête pour l’affichage de géométrie de type TIN devient alors:

select st_astext(ST_TinToGeomCollection(
st_geomFromText
            (
            'TIN(
            ((0 0 0, 1 0 0, 0 1 0, 0 0 0)),
            ((0 0 0, 1 0 0, 0 0 1, 0 0 0)),
            ((0 0 0, 0 0 1, 0 1 0, 0 0 0)),
            ((0 0 1, 1 0 0, 0 1 0, 0 0 1))
            )'

            ) ,0))
le second paramètre de la fonction ST_TinToGeomCollection étant le SRID de la géométrie, 4326 par exemple si vous êtes dans un système de coordonnées WGS84. Si vous ne le connaissez pas saisissez la valeur 0.

Affichage de données dans une base 3DCityDB contenant des données au format CityGML

Dans le panneau Query copiez la requête suivante :

WITH wallsurface as
(
SELECT a.id,c.root_id,st_collect(c.geometry) AS g,'255:128:0'::text AS rgb
FROM building a,thematic_surface b, surface_geometry c
WHERE b.building_id=a.id
AND b.type='WallSurface'
AND c.root_id=b.lod2_multi_surface_id
GROUP BY a.id,c.root_id
), roofsurface as
(
SELECT a.id,c.root_id,st_collect(c.geometry) AS g,'255:0:0'::text AS rgb
FROM building a,thematic_surface b, surface_geometry c
WHERE b.building_id=a.id
AND b.type='RoofSurface'
AND c.root_id=b.lod2_multi_surface_id
GROUP BY a.id,c.root_id
)
select id,root_id,g as geom,rgb
from wallsurface
union all
select id,root_id,g as geom,rgb
FROM roofsurface
ORDER BY id

image

Cliquez sur l’icone ‘Execute query

image
Dans la requête il est possible d’affecter une couleur à un type de géométrie, la syntaxe est le code R:G:B entre quotte suivi de AS rgb :
’255:128:0’ AS rgb  (wallsurface)
‘255:0:0’     AS rgb  (rootsurface)
image

 

Interface

   
image  
image Sortie de l’application.
image Connexion à une base de données.
   
image Panneau Requête
image Remet à blanc la fenêtre de requêtes
image Ouvre un fichier requête, qui ne doit contenir qu’une seule requête.
image Exécution de la requête.
image Collapse le panneau sur la gauche.
   
image Panneau affichage 3D
image Change la cuouleur du fond, Blanc ou Noir.
image Zoom Full.
image Zoom Plus.
image Zoom Moins.
image Efface le contenu de la fenêtre de rendu, toutes les géométries précédemment chargées sont supprimées.
   
image Panneau données attributaires
image Collapse le panneau sur le bas.
   
   
    

image onclusion.

Ce billet met à votre disposition un outil permettant de visualiser des géométries 3D stockées dans une base de données PostGIS 2.0 ou supérieure. Il me reste à implémenter l’affichage des points et des lignes.

Toutes suggestions d’amélioration seront les bienvenues, j’attends vos retours.

by Jérôme ROLLAND (noreply@blogger.com) at July 20, 2014 04:16 PM

A GeoSpatial World

PostGIS 3D Viewer introduction

 

image_thumb178 n viewer 3D pour PostGIS

image_thumb7[1] pg3DViewer_title6_thumb8
 

image_thumb175 ntroduction

Avec la version 2.0 PostGIS est entré dans le monde de la 3D, il est possible de stocker des géométries 3D avec l’ajout de la coordonnée Z dans toutes les géométries existantes et de stocker des géométries volumétriques avec l’introduction de deux nouveaux type  :

  • TIN(Triangulated Irregular Network, une collection de TRIANGLE)
image_thumb13
  • POLYHEDRALSURFACE (une collection de POLYGON)
image_thumb101

Cela fait déjà un an que la création d’un viewer pour visualiser les géométries 3D de PostGIS me trottait dans la tête suite à un échange avec Claus Nagel et Felix Kunde de virtualcitySYSTEMS à propos de l’outil d’import du format CityGML que j’avais réalisé CityGML vers PostGIS 2.0 , ils m’avaient conseillé alors de me diriger vers le développement d’un viewer 3D. Lorsque récemment j’ai repris contact avec eux pour leur parler de mon développement en cours, ils ont accueilli la nouvelle avec enthousiasme et ont bénéficié de versions beta, ce qui m’a grandement aidé pour mes choix dans la finalisation de mon application.

Aujourd’hui la 3D est de plus en plus présente dans le monde SIG avec le format CityGML et dans le monde de l’architecture avec le format IFC, importer ces formats dans une base de données PosgtgreSQL/PostGIS est une démarche naturelle car cela apporte toute la puissance de traitement des bases de données spatiale pour la gestion de ces données.

Pour le format CityGML la société virtualcitySYSTEMS met à disposition un schéma de base de données Open Source et gratuit pour le stockage et la gestion de ‘3D City Models’. Ce schéma de base de données est le reflet du modèle de données CityGML. Les données peuvent être importées via un outil d’import/export .

Affichage de données CityGML dans le viewer 3D

image_thumb14[6]

 

Continuer la lecture/Continue reading “PostGIS 3D Viewer”

by Jérôme ROLLAND (noreply@blogger.com) at July 20, 2014 04:12 PM

July 15, 2014

Stephen Mather

LiDAR and pointcloud extension pt 6

There’s all sorts of reasons why what I’ve been doing so far is wrong (see the previous 5 posts).  More on that later. But in case you needed 3D visualization of the chipping I’ve been working through, look no further:

Isometric 3D visualization of 3D chipping


by smathermather at July 15, 2014 12:32 AM

July 13, 2014

Boston GIS (Regina Obe, Leo Hsu)

PostGIS Minimalist X3D Viewer for PHP and ASP.NET

I've been thinking for a long time about creating an HTML based X3D viewer to visualize PostGIS 3D geometries. Now that I've started to document the new SFCGAL features, the need became more pressing. So I've created a very rudimentary one. You can check out the code on my github page https://github.com/robe2/postgis_x3d_viewer. I'm hoping to expand it in the future to allow people to pick materials for their geometries and also allow people to save the generated scene as a full X3D document.

x3d viewer for postgis

It utilizes the X3DOM Javascript library and JQuery for injecting X3D objects into the scene. I'm really excited about X3DOM JS. It was fairly trivial to work with and integrate the ST_AsX3D function. Utilizing it will help me stress test the ST_AsX3D function and patch up the holes in it. When I developed the function I was only interested in 3D geometries and didn't think too much about collections either. So as a result the function is kinda broken when working with 2D or Geometry Collections.

Features and Limitations

  • It has one scene and every query you do writes into the scene. This allows you to dump a bunch of objects on the same white board using different queries.
  • Similar to Raster Viewer I created a while ago, it's got a query log and a color picker. The color picker allows you to color each query output differently. I hope to augment this to allow picking different materials too.
  • Requires either .NET or PHP and of course a PostGIS 2+ enabled database
  • Similar to the Web raster viewer I created, it can only handle a query that outputs a single geometry (geometry mode) or an X3D scene snippet (raw mode)

Continue reading "PostGIS Minimalist X3D Viewer for PHP and ASP.NET"

by Regina Obe (nospam@example.com) at July 13, 2014 02:15 AM

July 12, 2014

Stephen Mather

LiDAR and pointcloud extension pt 5

Now for the crazy stuff:

Plaid-like overlay of vertical patches

The objective is to allow us to do vertical and horizontal summaries of our data. To do this, we’ll take chipped LiDAR input and further chip it vertically by classifying it.

First a classifier for height that we’ll use to do vertical splits on our point cloud chips:

CREATE OR REPLACE FUNCTION height_classify(numeric)
  RETURNS integer AS
$BODY$
 
SELECT
	CASE WHEN $1 < 10 THEN 1
	     WHEN $1 >= 10 AND $1 < 20 THEN 2
	     WHEN $1 >= 20 AND $1 < 30 THEN 3
	     WHEN $1 >= 30 AND $1 < 40 THEN 4
	     WHEN $1 >= 40 AND $1 < 50 THEN 5
	     WHEN $1 >= 50 AND $1 < 60 THEN 6
	     WHEN $1 >= 60 AND $1 < 70 THEN 7
	     WHEN $1 >= 70 AND $1 < 80 THEN 8
	     WHEN $1 >= 80 AND $1 < 90 THEN 9
	     WHEN $1 >= 90 AND $1 < 100 THEN 10
	     WHEN $1 >= 100 AND $1 < 110 THEN 11
	     WHEN $1 >= 110 AND $1 < 120 THEN 12
	     WHEN $1 >= 120 AND $1 < 130 THEN 13
	     WHEN $1 >= 130 AND $1 < 140 THEN 14
	     WHEN $1 >= 140 AND $1 < 150 THEN 15
	     WHEN $1 >= 150 AND $1 < 160 THEN 16
	     WHEN $1 >= 160 AND $1 < 170 THEN 17	     
	     WHEN $1 >= 170 AND $1 < 180 THEN 18
	     ELSE 0
	END
	FROM test;
 
$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

And now, let’s pull apart our point cloud, calculate heights from approximate ground, then put it all back together in chips grouped by height classes and aggregates of our original chips. Yes, I will explain more later. And don’t do this at home– I’m not even sure if it makes sense yet… :

/*
First we explode the patches into points,
 and along the way grab the minimum value
for the patch, and the id
*/
WITH lraw AS (
    SELECT PC_Explode(pa) AS ex,  PC_PatchMin(pa, 'Z') AS min, id
    FROM test1
    ),
/*
Calculate our height value ( Z - min )
*/
heights AS (
    SELECT PC_Get(ex, 'X') || ',' || PC_Get(ex, 'Y') || ',' -- grab X and Y
    || PC_Get(ex, 'Z') - min                -- calculate height
    || ',' ||
-- And return our remaining dimensions
    PC_Get(ex, 'Intensity') || ',' || PC_Get(ex, 'ReturnNumber') || ',' || PC_Get(ex, 'NumberOfReturns') || ',' ||
    PC_Get(ex, 'ScanDirectionFlag') || ',' || PC_Get(ex, 'EdgeOfFlightLine') || ',' ||
    PC_Get(ex, 'Classification') || ',' || PC_Get(ex, 'ScanAngleRank') || ',' || PC_Get(ex, 'UserData') || ',' ||
    PC_Get(ex, 'PointSourceId') || ',' || PC_Get(ex, 'Time') || ',' || PC_Get(ex, 'PointID') || ',' || PC_Get(ex, 'BlockID')
    AS xyheight, PC_Get(ex, 'Z') - min AS height, id FROM lraw
    ),
-- Now we can turn this aggregated text into an array and then a set of points
heightcloud AS (
    SELECT PC_MakePoint(1, string_to_array(xyheight, ',')::float8[]) pt, height_classify(height) heightclass, id FROM heights
    )
-- Finally, we bin the points back into patches grouped by our heigh class and original ids.
SELECT PC_Patch(pt), id / 20 AS id, heightclass FROM heightcloud GROUP BY id / 20, heightclass;

The resultant output should allow us to query our database by height above ground and patch. Now we can generate vertical summaries of our point cloud. Clever or dumb? It’s too early to tell.


by smathermather at July 12, 2014 01:56 AM

July 11, 2014

Stephen Mather

LiDAR and pointcloud extension pt 4

Height cloud overlayed with patchs

Now, we navigate into unknown (and perhaps silly) territory. Now we do point level calculations of heights, and drop the new calculated points back into the point cloud as though it were a Z value. I don’t recommend this code as a practice– this is as much as me thinking aloud as anything Caveat emptor:

/*
First we explode the patches into points,
 and along the way grab the minimum value
for the patch, and the id
*/
WITH lraw AS (
	SELECT PC_Explode(pa) AS ex,  PC_PatchMin(pa, 'Z') AS min, id
	FROM test1
	),
/*
Calculate our height value ( Z - min ) 
*/
heights AS (
	SELECT PC_Get(ex, 'X') || ',' || PC_Get(ex, 'Y') || ',' -- grab X and Y
	|| PC_Get(ex, 'Z') - min				-- calculate height
	|| ',' ||
-- And return our remaining dimensions
	PC_Get(ex, 'Intensity') || ',' || PC_Get(ex, 'ReturnNumber') || ',' || PC_Get(ex, 'NumberOfReturns') || ',' ||
	PC_Get(ex, 'ScanDirectionFlag') || ',' || PC_Get(ex, 'EdgeOfFlightLine') || ',' ||
	PC_Get(ex, 'Classification') || ',' || PC_Get(ex, 'ScanAngleRank') || ',' || PC_Get(ex, 'UserData') || ',' ||
	PC_Get(ex, 'PointSourceId') || ',' || PC_Get(ex, 'Time') || ',' || PC_Get(ex, 'PointID') || ',' || PC_Get(ex, 'BlockID')
	AS xyheight, id FROM lraw
	),
-- Now we can turn this aggregated text into an array and then a set of points
heightcloud AS (
	SELECT PC_MakePoint(1, string_to_array(xyheight, ',')::float8[]) pt, id FROM heights
	)
-- Finally, we bin the points back into patches
SELECT PC_Patch(pt) FROM heightcloud GROUP BY id / 20;

by smathermather at July 11, 2014 03:11 AM

July 09, 2014

Stephen Mather

LiDAR and pointcloud extension pt 3

Digging a little deeper. Ran the chipper on a smaller number of points and then am doing a little hacking to get height per chip (if you start to get lost, go to Paul Ramsey’s tutorial).

Here’s my pipeline file. Note the small chipper size– 20 points per chip.

<?xml version="1.0" encoding="utf-8"?>
<Pipeline version="1.0">
  <Writer type="drivers.pgpointcloud.writer">
    <Option name="connection">dbname='pggis' user='pggis' host='localhost'</Option>
    <Option name="table">test1</Option>
    <Option name="srid">3362</Option>
    <Filter type="filters.chipper">
      <Option name="capacity">20</Option>
      <Filter type="filters.cache">
        <Reader type="drivers.las.reader">
          <Option name="filename">60002250PAN.las</Option>
          <Option name="spatialreference">EPSG:3362</Option>
        </Reader>
      </Filter>
    </Filter>
  </Writer>
</Pipeline>

Easy enough to load (though slow for the sake of the chip size):

pdal pipeline las2pg.xml

Now we can do some really quick and dirty height calculations per chip — like what’s the difference between the minimum and maximum value in the chip:

DROP TABLE IF EXISTS test1_patches;

CREATE TABLE patch_heights AS
SELECT
  pa::geometry(Polygon, 3362) AS geom,
  PC_PatchAvg(pa, 'Z') AS elevation,
  PC_PatchMax(pa, 'Z') - PC_PatchMin(pa, 'Z') AS height
FROM test1;

Patch heights showing vegetation and open areas


by smathermather at July 09, 2014 02:55 AM

Stephen Mather

LiDAR and pointcloud extension pt 2

As much for my edification as anyones:
Want to get at just a few points within the pointcloud extension? A with query will get you there:

DROP TABLE IF EXISTS points;

CREATE TABLE points AS
WITH
-- get me some patches (our lidar points are stored in patches)
patches AS (
  SELECT pa FROM test
	LIMIT 2
),
-- Now, I want all the points from each patch
pa_pts AS (
  SELECT PC_Explode(pa) AS pts FROM patches
)
-- OK. Can I have that as standard old points, PostGIS style?
SELECT pts::geometry FROM pa_pts;
--k thanx

LiDAR points shown within their patches


by smathermather at July 09, 2014 02:35 AM

Stephen Mather

LiDAR and pointcloud extension

Paul Ramsey has a great tutorial on using the pointcloud extension with PostgreSQL / PostGIS: http://workshops.boundlessgeo.com/tutorial-lidar/

Image of lidar chipped into 400 point envelopes

You can get point cloud and all that goodness running a variety of ways. Probably the easiest is to download OpenGeo Suite from Boundless: http://boundlessgeo.com/solutions/opengeo-suite/download/

If you are an Ubuntu user, try a docker instance to run PostGIS with PDAL, pointcloud, etc in a container:

https://github.com/vpicavet/docker-pggis

Also, I’m working toward a simple installer script for Ubuntu 14.04 and for Vagrant, which is a fork of Vincent’s docker install:

https://github.com/smathermather/labr-pdal

One tip for those using, say QGIS instead of GeoServer for connecting and visualizing your data. Change any views created on the fly to tables, e.g.:

CREATE VIEW medford_patches AS
SELECT
  pa::geometry(Polygon, 4326) AS geom,
  PC_PatchAvg(pa, 'Z') AS elevation
FROM medford;

change to:

CREATE TABLE medford_patches AS
SELECT
pa::geometry(Polygon, 4326) AS geom,
PC_PatchAvg(pa, ‘Z’) AS elevation
FROM medford;

CREATE VIEW medford_patches AS
SELECT
  pa::geometry(Polygon, 4326) AS geom,
  PC_PatchAvg(pa, 'Z') AS elevation,
  id
FROM medford;

QGIS doesn’t like those views with casts for some reason… . without an id.

If the code I’m working on now works, it will do some height calculations using pointcloud, will show up on Github, and may include a little bonus pointcloud tutorial here. Stay tuned.


by smathermather at July 09, 2014 02:09 AM

July 07, 2014

Paul Ramsey

FOSS4G 2014 in Portland, Oregon, September 2014

Just a quick public service announcement for blog followers in the Pacific Northwest and environs: you've got a once in a not-quite-lifetime opportunity to attend the "Free and Open Source Software for Geospatial" (aka FOSS4G) conference this year in nearby Portland, Oregon, a city so hip they have trouble seeing over their pelvis.

Anyone in the GIS / mapping world should take the opportunity to go, to learn about what technology the open source world has available for you, to meet the folks writing the software, and the learn from other folks like you who are building cool things.

September 8th-13th, be there and be square.

by Paul Ramsey (noreply@blogger.com) at July 07, 2014 09:07 PM

July 03, 2014

Oslandia

QGIS versioning plugin

We developped a tool to manage data history, branches, and to work offline with your PostGIS-stored data and QGIS. Read more to get the insight of QGIS Versioning plugin.


The QGIS plugin is available in QGIS plugin repository, and you can fork it on GitHub too !


Capturing data offline

Introduction

Even if the necessity of data versioning often arises, no standard solution exist for databases.

The GeoGit project proposes a solution to store versioned geospatial data. There is also an existing plugin for QGIS, pgversion, which uses views and triggers to version a PostGIS database.

Unfortunately those solutions were not adapted to the specific constrains of this project, namely: using a PostGIS database as the main repository (excludes GeoGit) and the ability to working off-line (excludes pgversion).

The project we developed QGIS/PostGIS versioning looks like the following.



Design

The database is stored in a PostGIS schema, the complete schema is versioned (i.e. not individual tables). Revisions are identified by a revision number. A revision table in the versioned schema, called 'revisions', keeps track of the date, author, commit message and branch of all revisions.


Once a table structure is defined, three operations can be performed on rows: INSERT, DELETE and UPDATE. To be able to track history, every row is kept in the tables. Deleted rows are marked as such and updated rows are a combined insertion-deletion where the deleted and added rows are linked to one another as parent and child.


A total of five columns are needed for versioning the first branch:


PRIMARY KEY
a unique identifier across the table
branch_rev_begin
revision when this record was inserted
branch_rev_end
last revision for which this record exist (i.e. revision when it was deleted minus one)
branch_parent
in case the row has been inserted as the result of an update, this fields stores the hid of the row that has been updated
branch_child
in case the row has been marked as deleted as the result of an update, this field stores the hid of the row that has been inserted in its place.

For each additional branch, four additional columns are needed (the ones with the prefix branch_).


Note:
If the branch_rev_begin is null, it means that a row belongs to another branch.

SQL views are used to see the database for a given revision number. If we note 'rev' the revision we want to see. For each table, the condition for a row to be present is the view is:

(branch_rev_end IS NULL OR branch_rev_end >= rev) AND branch_rev_begin <= rev

In the special case of the current revision, or head revision, the condition reads:

branch_rev_end IS NULL AND branch_rev_begin IS NOT NULL
Note:
Since elements are not deleted (but merely marked as such) from an historized table, care must be taken with the definition of constrains, in particular the conceptual unicity of a field values.

Withing the PostGIS database, the views on revisions must be read-only and historized tables should not be edited directly. This is a basic principle for version control: editions must be made to working copies an then committed to the database. Please note that by default PostGIS 9.3 creates updatable views.


Workflow schema

This setup allows for multiple users to use and edit data offline from a central repository, and commit their modifications concurrently.



Working copies

Two kinds of working copies are available:

SpatiaLite working copies
They are meant to be used off-line. They consist of the versioned tables of a given versioned database (i.e. PostGIS schema) or any subset. For each table, only the elements that have not been marked as deleted in the head revision need to be present. Furthermore only a subset of the elements the user needs to edit can be selected (e.g. a spatial extend). To create a working copy (i.e. to checkout), tables from the versioned schema (or the aforementioned subsets) are converted to a SpatiaLite database using ogr2ogr.

PostGIS working copies
They are meant to be used when the connection to the original database will remain available. They are quite similar to pgversion working copies since they only store differences from a given revision (the one checked out).

The following description is aimed at understanding the inner workings of the qgis versioning plugin. The user does not need to perform the described operations manually.


For each versioned table in the working copy, a view is created with the suffix _view (e.g. mytable_view). Those views typically filters out the historization columns and shows the head revision. A set of triggers is defined to allow updating on those views (DELETE, UPDATE and INSERT).


The DELETE trigger simply marks the end revision of a given record.


The INSERT trigger create a new record and fills the branch_rev_begin field.


The UPDATE trigger create a new record and fills the branch_rev_begin and branch_parent fields. It then marks the parent record as deleted, and fills the branch_rev_end and branch_child fields.

Updating the working copy

Changes can be made to the database while editing the working copy. In order to reconcile those edition, the user needs to update the working copy.


When updating, a set of records can be in conflicts: the records for which the end revision has been set since the initial checkout or last update if any.


Multiple editions can be made to the same record. Therefore the child relation must be followed to the last child in order to present tu user with the latest state of a given conflicting feature.


Conflicts are stored in a table and identified with a conflict id and the tag 'theirs' or 'mine'. A DELETE trigger on this table is used for conflict resolution. On deletion of 'mine', the working copy edition is discarded, on deletion of 'theirs' the working copy edition is appended to the feature history (i.e. the working copy feature becomes a child of the last state of the feature in the historized database).

Committing the editions to the versionned database

If a working copy is up to date, the editions can be integrated in the versioned database. This operation consists simply in the insertion of a record in the revisions table, and, for each versioned table, the update of rows that are different and inserting rows that are not present.

Branching

A branch can be created from any revision by adding the four history columns and setting the branch_rev_begin field of features that are present in their revision.

Plugin interface tutorial

Groups are used for all versioning operations in QGIS since the versions are for a complete PostGIS schema or SpatiaLite database. The versioning toolbar will change depending on the group selected in the QGIS legend.

Note:
The group elements must share the same connection information (i.e. share the same database and schema for PostGIS working copies and revision views or share same SpatiaLite database for SpatiaLite working copies).

Versioning a PostGIS schema

Starting with an unversioned database, import a number of layers from a schema that needs to be versioned into a QGIS project. Once the layers are imported, they must be grouped together.



Selecting the newly created group will cause the versioning toolbar to display the historize button (green V). On click a confirmation is requested to version the database schema.



The versioned layers are imported in a new group and the original layers are removed from the project.


Note:
The symobology is not kept in the process.


Working with a versioned PostGIS schema

Versioned layers can be imported in QGIS. The layers must be from a head revision or a view on any revision.



Once the layers are in QGIS, they must be grouped.



For PostGIS groups at head revision, the versioning plugin allows the user to create a SpatiaLite or a PostGIS working copy, create a view on a given revision or create a branch. A corresponding group will be imported in QGIS.



If the user chooses to create a SpatiaLite working copy, he will be asked to select a file to store the working copy.



Commiting changes

Once the working copy is imported in QGIS, the user can start edition of its layers. With SpatiaLite working copies, this edition can be done off-line.



When the user is done with edition, he can commit the changes to the database and if commit is feasible (i.e. the working copy is up to date with the versioned database), he will be prompted for a commit message and subsequently be informed of the revision number he committed.




If the commit is not feasible, the user will be informed that he must update his working copy prior to commit.


Resolving conflicts

Conflicts are detected during update, the user is informed, and conflicts layers are imported into QGIS.



To resolve conflicts, the user can open the conflict layer's attribute table. Selected entries are also selected on the map canvas and the user can decide which version, either his or the database's, he wants to keep. User version is tagged with 'mine' and database version with 'theirs'. The conflict is resolved by deleting the unwanted entry in the conflict layer.


Note:
On deletion of one conflict entry, both entries are removed (by a trigger) but the attribute table (and canvas) are not refreshed. As a workaround, the user can close and re-open the attribute table to see the actual state of the conflict table.

Once the conflict table is empty, the commit can be done.



Restrictions

Due to design choices and tools used for conversion to SpatiaLite, a number of restrictions apply to the versioned database:


  • schemas, tables and branch names should not have space, caps or quotes
  • tables must have primary keys
  • columns are lowercase (because of conversion to SpatiaLite) but can have spaces (not that it's recommended)
  • geometry column is geom in PostGIS, GEOMETRY in SpatiaLite

Note
Do not edit OGC_FID or ROWID
Note
The constrains on the tables are be lost in the PostGIS to SpatiaLite conversion.

Known bug

The conflict layer won't be loaded automatically is it has no geometry. The user will have to load it manually.

by Vincent Mora at July 03, 2014 10:00 AM

June 30, 2014

Oslandia

QGIS 2.4 release out, Oslandia inside

There is a new QGIS release out : version 2.4, codename Chugiak is now available. Binary packages for your platform have been generated, and you can directly download and try out this new release of the famous Desktop GIS software.



QGIS 2.4 has a lot of new features in all of its components. There is a visual changelog available where you can discover QGIS improvements.


Oslandia inside

Oslandia is a QGIS core contributor, and we have been busy improving QGIS 2.4. We contributed to various of these new features. Here are a few enhancements we developped.

Predefined scales mode for atlas maps

When working with atlas map items, you can now specify a predefined scale mode for the map. It will use the best fitting option from the list of predefined scales in you your project properties settings (see Project -> Project Properties -> General -> Project Scales to configure these predefined scales).



This feature has been funded by the city of Uster

New Inverted Polygon renderer

The biggest feature Oslandia developped is the inverted Polygon renderer. This feature has been funded by Agence de l'Eau Adour-Garonne and mainly developped by Hugo Mercier.


A new renderer has been added for polygon features, allowing you to style everything outside your polygons. This can be useful for highlighting areas, or for creating cartographic masks. When used with new shapeburst style, you can now produce output as shown in the image for this entry.


New Mask Plugin

Alongside with the inverted polygon renderer feature, Oslandia developped a new Mask plugin. It enables you to make atlases focusing on the specific feature you are interested in, occulting the rest with a really nice effect. Furthermore, it helps masking the labels when you mask geometry objects.


You can read more details on these features on this blog post.


This plugin has also be funded by Agence de l'Eau Adour Garonne.

Layered SVG export

Another feature implemented in this version is the ability to export layered SVG files. Beforehand, all features were exported as a single layer, whatever the QGIS layer was. Now you can use Inkscape or Illustrator and their layering capabilities to finish the design of your map with greater ease of use. There also is an option to vectorize labels.


This feature has been funded by Agence de Développement du Grand Amiénois (ADUGA) .

WAsP format support

The WAsP format is the standard format for roughness and elevation in the wind science field. This format was not supported by QGIS until recently, when Vincent Mora added WAsP to QGIS supported GIS file formats. In fact, we did better as we implemented WAsP support in GDAL/OGR , so that any software using this library is now able to read and write WASP files. WAsP is available starting from GDAL/OGR >= 1.11.0.


This was an opportunity to add Vincent Mora as an official GDAL/OGR commiter, in charge of maintaining this driver. This feature will enable wind management operations to be completed on QGIS with a better user experience. No more file conversion before working on the GIS side. We also developped a companion plugin to manage data simplification when needed. It is available in QGIS plugins repository.

With this work, QGIS becomes a great complement to opensource computational wind engineering softwares like ZephyTools .


This work has been funded by La Compagnie du Vent

Epanet Plugin

Oslandia has integrated the EPANET water distribution model inside QGIS Processing, as a plugin. You can read more on this plugin on this blog .


Epanet integration has been funded by European funds and the GIS office of Apavil, Romania.

Vizitown Plugin

Vizitown is part of our efforts on 3D GIS development , alongside PostGIS 3D and more. It is a QGIS plugin allowing users to display QGIS layers in 3D in a Three.js / WebGL environment, in a browser. It can leverage PostGIS 3D, and display live data from the database, as well as other sources of data. It can display DEM, a raster background, 2D vector data draped on the DEM, 2.5D data (e.g. buildings), or real 3D Meshes. The user can set a symbology in QGIS and see the modifications live in the browser in 3D.


You can see Vizitown in action on Youtube . Vizitown has been developped with IG3 students from ESIPE .

Multiple bugfixes

Oslandia also work continuously on improving QGIS quality, and we try to fix as many bugs as we can.


These bugfixes are funded by our QGIS support offer clients, and also by the french Ministère de l'environnement and Agence de l'Eau Adour-Garonne.

What next ?

We continue to work on new QGIS features, corrections, refactoring and integration with other tools. We namely studied a possible unification of all database-like features in QGIS, using SQLITE/Spatialite features. We intend to work on Native Read+Write support of Mapinfo TAB files.


We offer a wide range of services around QGIS, be it for training, assistance, development, or consulting in general.


We also propose various support opportunities for QGIS . This is the best way for you to improve the quality of this software, contribute to its development, and ensure that you can work in good conditions without having to worry about potential bugs. Our team of experienced engineers, who contribute to QGIS core, will be available in any case of critical bug.


We can offer you personalized support, with specific conditions and fares. Do not hesitate to contact us at infos@oslandia.com .

by Vincent Picavet at June 30, 2014 12:00 PM

June 25, 2014

shisaa.jp

Postgis, PostgreSQL's spatial partner - Part 3

You have arrived at the final chapter of this PostGIS introduction series. Before continuing, I recommend you read chapter one and chapter two first.

In the last chapter we finished by doing some real world distance measuring and we saw how different projections pushed forward different results.

Today I would like to take this practical approach a bit further and continue our work with real world data by showing you around the town of Kin in Okinawa. The town where I live.

A word before we start

In this chapter I want to do a few experiments together with you on real world data. To gather this data, I would like to use OpenStreetMap because it is not only open but also gives us handy tools to export map information.

We will use a tool called osm2pgsql to load our OSM data into PostGIS enable tables.

However, it is more common to import and export real world GIS data by using the semi-closed ESRI standard shapefile format. OpenStreetMap does not support exporting to this shapefile format directly, but exports to a more open XML file (.osm) instead.

Therefor, near the end of this post, we will briefly cover these shapefiles as well and see how we could import them into our PostgreSQL database. But for the majority of our work today, I will focus on the OpenStreetMap approach.

The preparation

Let us commence with this adventure by first getting all the GIS data related to the whole of Okinawa. We will only be interested in the data related to Kin town, but I need you to pull in a data set that is large enough (but still tiny in PostgreSQL terms) for us to experiment with indexing.

Hop online and download the file being served at the following URL: openstreetmap.org Okinawa island It is a file of roughly 180 Mb and covers most of the Okinawan main island. Save the presented "map" file.

Next we will need to install a third party tool which is specifically designed to import this OSM file into PostGIS. This tool is called osm2pgsql and is available in many Linux distributions.

On a Debian system:

apt-get install osm2pgsql

Loading foreign data

Now we are ready to load in this data. But first, let us clean our "gis" database we used before.

Since all these import tools will create their own PostGIS enabled tables, we can delete our "shapes" table. Connect to your "gis" database and drop this table:

DROP TABLE shapes;

Using this new tool, repopulate the "gis" database with the data you just downloaded:

osm2pgsql -s -U postgres -d gis map

If everything went okay, you will get a small report containing the information about all the tables and indexes that where created.

Let us see what we just did.

First we ran osm2pgsql with the -s flag. This flag enabled slim mode, which means it will use a database on disk, rather then processing all the GIS data in RAM. The latter does not only potentially slow down your machine for larger data sets, but it enables less features to be available.

Next we tell the tool to connect as the user "postgres" and load the data into the "gis" database. The final argument is the "map" file you just downloaded.

What do we have now?

Open up a database console and let us describe our database to see what this tool just did:

\d

As you can see, it inserted 7 new tables:

Schema |        Name        | Type  |  Owner   
--------+--------------------+-------+----------
public | geography_columns  | view  | postgres
public | geometry_columns   | view  | postgres
public | planet_osm_line    | table | postgres
public | planet_osm_nodes   | table | postgres
public | planet_osm_point   | table | postgres
public | planet_osm_polygon | table | postgres
public | planet_osm_rels    | table | postgres
public | planet_osm_roads   | table | postgres
public | planet_osm_ways    | table | postgres
public | raster_columns     | view  | postgres
public | raster_overviews   | view  | postgres
public | spatial_ref_sys    | table | postgres

The other 5 views and tables are the good old PostGIS bookkeeping.

It is also important, yet less relevant for our work here today, to know that these tables, or rather the way osm2pgsql imports, is optimized to work with Mapnik. Mapnik is an open-source map rendering software package used for both web and offline usage.

The tables that are imported contain many different types of information. Let me quickly go over them to give you a basic feeling of how the import happened:

  • planet_osm_line: holds all non-closed pieces of geometry (called ways) at a high resolution. They mostly represent actual roads and are used when looking at a small, zoomed-in detail of a map.
  • planet_osm_nodes: an intermediate table that holds the raw point data (points in lat/long) with a corresponding "osm_id" to map them to other tables
  • planet_osm_point: holds all points-of-interest together with their OSM tags - tags that describe what they represent
  • planet_osm_polygon: holds all closed piece of geometry (also called ways) like buildings, parks, lakes, areas, ...
  • planet_osm_rels: an intermediate table that holds extra connecting information about polygons
  • planet_osm_roads: holds lower resolution, non-closed piece of geometry in contrast with "planet_osm_line". This data is used when looking at a greater distance, covering much area and thus not much detail about smaller, local roads.
  • planet_osm_ways: an intermediate table which holds non-closed geometry in raw format

We will now continue working with a small subset of this data.

Let us take a peek at the Polygons tables for example. First, let us see what we have available:

\d planet_osm_polygon

That is quite a big list, but the major part of these columns are of mere TEXT type and contain human information about the geometry stored. These columns corresponds with the way OpenStreetMap categorizes their data and with the way you could use the Mapnik software described above.

Let us do a targeted query:

SELECT name, building, ST_AsText(way)
    FROM planet_osm_polygon
    WHERE building = 'industrial';

Notice that I use the output function ST_AsText() to convert to a human readable WKT string. Also, I am only interested in some of the industrial buildings, so I set the building type to industrial.

The result:

                    name                      |  building  |                                                                     st_astext                                                                   
----------------------------------------------+------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------
 沖ハム (Okiham)                               | industrial | POLYGON((14221927.83 3049797.01,14222009.77 3049839.68,14222074.84 3049714.68,14222028.9 3049690.76,14221996.33 3049753.33,14221960.34 3049734.58,14221927.83 3049797.01))
Kin Thermal Power Plant Coal storage building | industrial | POLYGON((14239931.42 3054117.72,14239990.49 3054224.25,14240230.15 3054091.38,14240171.08 3053984.84,14239931.42 3054117.72))
Kin Thermal Power Plant Exhaust tower         | industrial | POLYGON((14240167.1 3054497.14,14240172.26 3054507.93,14240176.04 3054515.82,14240195.76 3054506.39,14240186.84 3054487.7,14240167.1 3054497.14))

We get back three records containing one industrial building each, described with a closed polygon. Cool.

Now, I can assure you that Okinawa has more then three industrial buildings, but do remember that we are looking at a rather rural island. OpenStreetMap relies greatly on user generated content and there simply are not many users who have felt the need to index the industrial buildings here in this neck of the woods.

The planet_osm_polygon table does contain little over 6000 buildings of various types, which is still a small number, but for our purpose today I am only interested in the latter two, which both lie here in Kin town.

Also, if you would, for example, take a chunk of Tokyo, where there are hundreds of active OpenStreetMap contributors, you will find that many buildings are present and are sometimes even more accurately represented then some other online proprietary mapping solutions offered by some famous search engines. Ahum.

Before continuing, though, I would like to delete two GiST indexes that "osm2pgsql" made for us, purely to be able to demonstrate the importance of an index.

For now, just take my word and delete the indexes on all the geometry columns of the tables we will use today:

DROP INDEX planet_osm_line_index;
DROP INDEX planet_osm_polygon_index;

Then perform a VACUUM:

VACUUM ANALYZE planet_osm_line;
VACUUM ANALYZE planet_osm_polygon;

VACUUM together with ANALYZE will force PostgreSQL to recheck the whole table for any changed conditions, as is the case since we removed the index.

The first thing I would like to find out is how large these building actually are. We cannot measure how tall they are, for we are working with two dimensional data here, but we can measure their footprint on the map.

Since PostGIS makes all of our work easy, we could simply employ a function to tell us this information:

SELECT ST_Area(way)
    FROM planet_osm_polygon
    WHERE building = 'industrial';

And get back:

    st_area      
------------------
10155.3935499731
33381.1043500491
452.9464999972

As we know from the previous chapter, to be able to know what these numbers mean, we have to find out in which SRID this data was saved. You could either describe the table again and look at the geometry column description, or use an accessor function ST_SRID(), to find it:

SELECT ST_SRID(way)
    FROM planet_osm_polygon;
    WHERE building = 'industrial';

And get back:

 st_srid 
---------
900913
900913
900913

You could also query the PostGIS bookkeeping directly and look in the geometry_columns view:

SELECT f_tablename, f_geometry_column, coord_dimension, srid, type
    FROM geometry_columns;

This view holds information about all the geometry columns in our PostGIS enabled database. Our above query will return a list containing all the GIS describing information we saw in the previous chapter.

Nice. Both our buildings are stored in a geometry column and have an SRID of 900913. We can now use our spatial_ref_sys table to look up this ID:

SELECT srid, auth_name, auth_srid, srtext, proj4text
    FROM spatial_ref_sys
    WHERE srid = 900913;

As you can see, this is basically a Mercator projection used by OpenStreetMap. In the "proj4text" column we can see that its units are meters.

This thus means that the information we get back is in square Meters.

In this map (only looking at the latter two Kin buildings) we thus have a building with a total area of 33 square Kilometers and a more modest building of around 452 square Meters. The former is a coal storage facility belonging to the Kin Thermal Power Plant and is indeed huge. The second building represents the exhaust tower of that same plant.

You have just measured the area these buildings occupy, very neat right?

Now, let us find out which road runs next to this power plant, just in case we wish to drive to there. It is important to note that OSM (and many other mapping solutions) divide roads into different types.

You have trunk roads, highways, secondary roads, tertiary roads, etc. I am now interested to find the nearest secondary road.

To get a list of all the secondary roads in Okinawa, simply query the planet_osm_roads table:

SELECT ref, ST_AsText(way) 
    FROM planet_osm_line
    WHERE highway = 'secondary';

We now get back all the linestring objects together with their reference inside of OSM. The reference refers to the actual route number each road has.

The total count should be around 3215 pieces of geometry, which is already a nice list to work with.

Let us now see which of these roads is closest to our coal storage building.

To find out how far something is (nearest neighbor search) we could use our ST_Distance() function we used in the previous chapter and perform the following lookup:

SELECT road.highway, road.ref, ST_Distance(road.way, building.way) AS distance
    FROM planet_osm_polygon AS building, planet_osm_line AS road
    WHERE road.highway = 'secondary' AND building.name = 'Kin Thermal Power Plant Coal storage building'
    ORDER BY distance;

This will bring us:

 highway  | ref |     distance     
-----------+-----+------------------
secondary | 329 | 417.374986575458
secondary | 104 | 2258.90394593648
secondary | 104 | 2709.00178089638
secondary | 104 | 2745.76782385198
secondary | 234 | 5897.78205314507
...

Cool, secondary route 329 is the closest to our coal storage building with a distance of 417 meters.

While this will return quite accurate results, there is one problem with this query. Indexes are not being used. And every time an index is potentially left alone, you should start to worry, especially with larger data sets.

How do I know they are ignored? Simple, we did not make any indexes (and we deleted the ones made by "osm2pgsql")...which makes me pretty sure we cannot use them.

I refer you to chapter three of my PostgreSQL Full Text series where I talk a bit more about GiST and B-Tree index types. And, as I also say in that chapter, I highly recommend reading Markus Winand's Use The Index, Luke series, which explains in great detail how database indexes work.

The first thing to realize is that an index will only be used if the data set on which it is build is of sufficient size. PostgreSQL has an AI build in, called the query planner, which will make a decision on whether or not to use an index.

If your data set is small enough a more traditional Sequential Scan will be faster or equal.

To know what is going on exactly and to know how fast our query runs, we have the EXPLAIN command at our disposal.

Speeding things up

Let us EXPLAIN the query we have just run:

EXPLAIN ANALYZE SELECT road.highway, road.ref, ST_Distance(road.way, building.way) AS distance
    FROM planet_osm_polygon AS building, planet_osm_line AS road
    WHERE road.highway = 'secondary' AND building.name = 'Kin Thermal Power Plant Coal storage building'
    ORDER BY distance;

We simply put the keyword EXPLAIN (and ANALYZE to give us total runtime) right in front of our normal query.

The result:

Sort  (cost=5047.50..5055.32 rows=3129 width=391) (actual time=41.481..41.815 rows=3215 loops=1)
    Sort Key: (st_distance(road.way, building.way))
    Sort Method: quicksort  Memory: 348kB
    ->  Nested Loop  (cost=0.00..4309.34 rows=3129 width=391) (actual time=1.188..38.617 rows=3215 loops=1)
        ->  Seq Scan on planet_osm_polygon building  (cost=0.00..279.01 rows=1 width=207) (actual time=0.981..1.409 rows=1 loops=1)
            Filter: (name = 'Kin Thermal Power Plant Coal storage building'::text)
            Rows Removed by Filter: 6320
        ->  Seq Scan on planet_osm_line road  (cost=0.00..3216.79 rows=3129 width=184) (actual time=0.166..26.524 rows=3215 loops=1)
            Filter: (highway = 'secondary'::text)
            Rows Removed by Filter: 73488
Total runtime: 42.153 ms

That is a lot of output, but it shows you how the internal planner executes our query and which decisions it makes along the way.

To fully interpret a query plan (this is still a simple one), a lot more knowledge is needed and this would easily deserve its own series. I am by far not an expert in the query planner (though it is an interesting study topic), but I will do my best to extract the important bits we need for our direct performance tuning.

A query plan is always made up out of nested nodes, the parent node containing all the accumulated information (costs, rows, ...) of its child nodes.

Inside the nested loop parent node we see above, we can find that the planner decided to use two filters, which correspond to the WHERE clause conditions of our query (building.name and road.highway). You can see that both child nodes are of Seq Scan type, which means Sequential Scan. These types of nodes scan the whole table, simply from top to bottom, directly from disk.

Another important thing to note is the total time this query costs, which is 42.153 ms. The time reported here is the time on my local machine, depending on how decent your computer is, this time could vary.

A detail not to forget when looking at this timing, is the fact that it is slightly skewed if compared to real-world application use:

  • We neglect network/client traffic. This query now runs internally and does not need to communicate with a client driver (which almost always brings extra overhead)
  • The time measurement itself also introduces overhead.

The total runtime from our above plan does not sound as a big number, but we are working with a rather small data set - the area of Okinawa is large, but the geometry is rather sparse.

So our first reaction should be: this can be better.

First, let us try to get rid of these sequential scans, for they are a clear indication that the planner does not use an index.

Creating indexes

In our case we want to make two types of indexes:

  • Indexes on our "meta" data, the names and other attributes describing out geometrical data
  • Indexes that actually index our geometrical data itself

Let us start with our attributes columns.

These are all simple VARCHAR, TEXT or INT columns, so the good old Balanced Tree or B-Tree can be used here. In our query above we use "road.highway" and "building.name" in our lookup, so let us make a couple of indexes that adhere to this query. Remember, an index only makes sense if it is built the same way your queries question your data.

First, the "highway" column of the "planet_osm_line" table:

CREATE INDEX planet_osm_line_highway_index ON planet_osm_line(highway);

The syntax is trivial. You simply tell PostgreSQL to create an index, give it a name, and tell it on which column(s) of which table you want it to be built. PostgreSQL will always default to the B-Tree index type.

Next, the name column:

CREATE INDEX planet_osm_polygon_name_index ON planet_osm_polygon(name);

Now perform another VACUUM ANALYZE on both tables:

VACUUM ANALYZE planet_osm_line;
VACUUM ANALYZE planet_osm_polygon;

Let us run explain again on the exact same query:

Sort  (cost=4058.73..4066.56 rows=3129 width=394) (actual time=20.817..21.149 rows=3215 loops=1)
    Sort Key: (st_distance(road.way, building.way))
    Sort Method: quicksort  Memory: 348kB
    ->  Nested Loop  (cost=72.95..3310.07 rows=3129 width=394) (actual time=1.356..17.743 rows=3215 loops=1)
        ->  Index Scan using planet_osm_polygon_name_index on planet_osm_polygon building  (cost=0.28..8.30 rows=1 width=207) (actual time=0.054..0.056 rows=1 loops=1)
            Index Cond: (name = 'Kin Thermal Power Plant Coal storage building'::text)
        ->  Bitmap Heap Scan on planet_osm_line road  (cost=72.67..2488.23 rows=3129 width=187) (actual time=1.258..4.661 rows=3215 loops=1)
            Recheck Cond: (highway = 'secondary'::text)
            ->  Bitmap Index Scan on planet_osm_line_highway_index  (cost=0.00..71.89 rows=3129 width=0) (actual time=0.864..0.864 rows=3215 loops=1)
                   Index Cond: (highway = 'secondary'::text)
Total runtime: 21.527 ms

You can see that we now traded our Seq Scan for Index Scan and Bitmap Heap Scan, which indicates that our attribute indexes are being used, yay!

The so-called Bitmap Heap Scan, instead of a Sequential Scan, is performed when the planner decides it can use the index to gather all the rows it thinks it needs, sort them in logical order and then fetch the data from the table on disk in the most optimized way possible (trying to open each disk page only once).

The order by which the Bitmap Heap Scan arranges the data is directed by the child node aka the Bitmap Index Scan. This latter type of node is the one doing the actual searching inside the index. Because in our WHERE clause we have a condition which tells PostgreSQL to limit the rows to the ones of "highway" type "secondary", the Bitmap Index Scan fetches the needed rows from our B-Tree index we just made and passes them to its parent, the Bitmap Heap Scan, which then goes on to order the geometry rows to be fetched.

This already helped much, for our query runtime dropped to half. Now, let us make the indexes for our actual geometry, and see the effect:

CREATE INDEX planet_osm_line_way ON planet_osm_line USING gist(way);
CREATE INDEX planet_osm_polygon_way ON planet_osm_polygon USING gist(way);

Creating a GiST index is quite similar to a normal B-Tree index. The only difference here is that you specify the index to be build with GiST.

Vacuum:

VACUUM ANALYZE planet_osm_line;
VACUUM ANALYZE planet_osm_polygon;

Now poke it again with the same query and see our new plan:

Sort  (cost=4038.82..4046.54 rows=3089 width=395) (actual time=21.137..21.479 rows=3215 loops=1)
    Sort Key: (st_distance(road.way, building.way))
    Sort Method: quicksort  Memory: 348kB
    ->  Nested Loop  (cost=72.64..3299.76 rows=3089 width=395) (actual time=1.382..17.858 rows=3215 loops=1)
        ->  Index Scan using planet_osm_polygon_name_index on planet_osm_polygon building  (cost=0.28..8.30 rows=1 width=207) (actual time=0.041..0.044 rows=1 loops=1)
            Index Cond: (name = 'Kin Thermal Power Plant Coal storage building'::text)
        ->  Bitmap Heap Scan on planet_osm_line road  (cost=72.36..2488.32 rows=3089 width=188) (actual time=1.297..4.726 rows=3215 loops=1)
            Recheck Cond: (highway = 'secondary'::text)
            ->  Bitmap Index Scan on planet_osm_line_highway_index  (cost=0.00..71.59 rows=3089 width=0) (actual time=0.866..0.866 rows=3215 loops=1)
                  Index Cond: (highway = 'secondary'::text)
Total runtime: 21.873 ms

Hmm, the plan did not change at all, and our runtime is roughly identical. Why is our performance still the same?

The culprit here is ST_Distance().

As it turns out, this function is unable to use the GiST index and is therefor not a good candidate to set loose on your whole result set. The same goes for the ST_Area() function, by the way.

So we need a way to limit the amount of records we do this expensive calculation on.

ST_DWithin()

We introduce a new function: ST_DWithin(). This function could be our savior in this case, for it does use the GiST index.

Whether or not a function (or operator) can use the GiST index, depends on if it uses bounding boxes when performing calculations. The reason why is because GiST indexes mainly store bounding box information and not the exact geometry itself.

ST_DWithin() checks if given geometry is within a radius of another piece of geometry and simply returns TRUE or FALSE. We can thus use it in our WHERE clause to filter out geometry for which it returns FALSE (and thus not falls within the radius). It performs this check using bounding boxes, and thus is able to retrieve this information from our GiST index.

Let me present you with a query that limits the result set based on what ST_DWithin() finds:

EXPLAIN ANALYZE SELECT road.highway, road.ref, ST_Distance(road.way, building.way) AS distance
    FROM planet_osm_polygon AS building, planet_osm_line AS road
    WHERE road.highway = 'secondary'
        AND building.name = 'Kin Thermal Power Plant Coal storage building'
        AND ST_DWithin(road.way, building.way, 10000)
    ORDER BY distance;

As you can see we simply added one more WHERE clause to limit the returned geometry by radius. This will result in the following plan:

Sort  (cost=45.66..45.67 rows=1 width=395) (actual time=6.048..6.052 rows=27 loops=1)
    Sort Key: (st_distance(road.way, building.way))
    Sort Method: quicksort  Memory: 27kB
    ->  Nested Loop  (cost=4.63..45.65 rows=1 width=395) (actual time=3.157..6.005 rows=27 loops=1)
        ->  Index Scan using planet_osm_polygon_name_index on planet_osm_polygon building  (cost=0.28..8.30 rows=1 width=207) (actual time=0.051..0.054 rows=1 loops=1)
            Index Cond: (name = 'Kin Thermal Power Plant Coal storage building'::text)
        ->  Bitmap Heap Scan on planet_osm_line road  (cost=4.34..37.09 rows=1 width=188) (actual time=3.090..5.771 rows=27 loops=1)
            Recheck Cond: (way && st_expand(building.way, 10000::double precision))
            Filter: ((highway = 'secondary'::text) AND (building.way && st_expand(way, 10000::double precision)) AND _st_dwithin(way, building.way, 10000::double precision))
            Rows Removed by Filter: 4838
                ->  Bitmap Index Scan on planet_osm_line_way  (cost=0.00..4.34 rows=8 width=0) (actual time=1.978..1.978 rows=4865 loops=1)
                    Index Cond: (way && st_expand(building.way, 10000::double precision))
Total runtime: 6.181 ms

Good. We have just gone down to only 6.181 ms, That seems to be much more efficient.

As you can see, our query plan got a few new rows. The main thing to notice is the fact that our Bitmap Heap Scan got another Recheck Cond, our expanded ST_DWithin() condition. More to the bottom, you can see that the condition is being pulled from the GiST index:

Index Cond: (way && st_expand(building.way, 10000::double precision))

This seems to be a much more desirable and scalable query.

But there is a drawback, though ST_DWithin() will make for speedy results, it works only by giving it a fixed radius.

As you can see from our usage, we call the function as follows: ST_DWithin(road.way, building.way, 10000). The last argument, "10000", tells us how big the search radius is. In this case our geometry is in meters, so this means we search in a radius of 10 Km.

This static radius number is quite arbitrary and might not always be desirable. What other options do we have without compromising performance too much?

Operators

Another addition of PostGIS we have not talked about much up until now are the spatial operators we have available. You have a total of 16 operators you can use to perform matches on your GIS data.

You have straightforward operators like &&, which returns TRUE if one piece of geometry intersects with another (bounding box calculation) or the << which returns TRUE if one object is fully to the left of another object.

But there are more interesting ones like the <-> and the <#> operators.

The first operator, <->, returns the distance between two points. If you feed it other types of geometry (like a linestring of polygon) it will first draw a bounding box around that geometry and perform a point calculation by using the bounding box centroids. A centroid is the calculated center of a piece of geometry (the drawn bounding box in our case).

The second, <#>, acts completely the same, but works directly on bounding boxes of given geometry. In our case, since we are not working with points, it would make more sense to use this operator.

The big advantage of this distance calculation operator is, once more, the fact that it too calculates using a bounding box and is thus able to use a GiST index. However, the ST_Distance() function calculates distances by finding two points on the given geometry most close to each other, which serves the most accurate result. The <#> operator, as said before, stretches a bounding box around each piece of geometry and therefor deforms our objects, making for less accurate distance measuring.

It is therefor not wise to use <#> to calculate accurate distances, but it is a life saver to sort away geometry that is too far away for our interest.

So a proper usage would be to first roughly limit the result set using the <#> operator and then more accurately measure the distance of, say, the first 50 matches with our famous ST_Distance().

Before we can continue, it is important to point out that both the <-> and <#> operator can only use the GiST index when either the left or right hand side of the operator is a constant or fixed piece of geometry. This means we have to provide actual geometry using a constructor function.

There are other ways around this limitation by, for example as Alexandre Neto points out on the PostGIS mailing list, providing your own function which converts our "dynamic" geometry into a constant.

But this would make this post run way past its initial focus. Let us simply try by providing a fixed piece of geometry. The fixed piece is, of course, still our "Kin Thermal Power Plant Coal storage building", but converted into WKT:

EXPLAIN ANALYZE WITH distance AS (
    SELECT way AS road, ref AS route
        FROM planet_osm_line
        WHERE highway = 'secondary'
        ORDER BY ST_GeomFromText('POLYGON((14239931.42 3054117.72,14239990.49 3054224.25,14240230.15 3054091.38,14240171.08 3053984.84,14239931.42 3054117.72))', 900913) <#> way
        LIMIT 50
        )
SELECT ST_Distance(ST_GeomFromText('POLYGON((14239931.42 3054117.72,14239990.49 3054224.25,14240230.15 3054091.38,14240171.08 3053984.84,14239931.42 3054117.72))', 900913), road) AS true_distance, route
    FROM distance
    ORDER BY true_distance
    LIMIT 1;

This query uses a Common Table Expression or CTE (you could also use a simpler subquery) to first get a rough result set of about 50 rows based on what <#> finds. Then only on those 50 rows do we perform our more expensive, index-agnostic distance calculation.

This results in the following plan and runtime:

Limit  (cost=274.57..274.57 rows=1 width=64) (actual time=11.236..11.237 rows=1 loops=1)
    CTE distance
        ->  Limit  (cost=0.28..260.82 rows=50 width=173) (actual time=0.389..10.764 rows=50 loops=1)
            ->  Index Scan using planet_osm_line_way on planet_osm_line  (cost=0.28..16362.19 rows=3140 width=173) (actual time=0.389..10.745 rows=50 loops=1)
                Order By: (way <#> '010300002031BF0D000100000005000000D7A3706D17296B41C3F528DC124D47417B14AECF1E296B4100000020484D4741CDCCCCC43C296B410AD7A3B0054D4741295C8F6235296B41B81E856BD04C4741D7A3706D17296B41C3F528DC124D4741'::geometry)
                Filter: (highway = 'secondary'::text)
                Rows Removed by Filter: 4562
            ->  Sort  (cost=13.75..13.88 rows=50 width=64) (actual time=11.234..11.234 rows=1 loops=1)
                Sort Key: (st_distance('010300002031BF0D000100000005000000D7A3706D17296B41C3F528DC124D47417B14AECF1E296B4100000020484D4741CDCCCCC43C296B410AD7A3B0054D4741295C8F6235296B41B81E856BD04C4741D7A3706D17296B41C3F528DC124D4741'::geometry, distance.road))
                Sort Method: top-N heapsort  Memory: 25kB
    ->  CTE Scan on distance  (cost=0.00..13.50 rows=50 width=64) (actual time=0.412..11.188 rows=50 loops=1)
Total runtime: 11.268 ms

As you can see, we are now using the GiST index "planet_osm_line_way", which was what we were after.

This yields roughly the same runtime as with our ST_DWithin(), but without the arbitrary distance setting. We indeed have a somewhat arbitrary limiter of 50, but this is much less severe then a distance limiter.

Even if the closest secondary road is 100 Km from our building, the above query would still find it whereas our previous query would return nothing.

One more for the road home

Let us do a few more fun calculations on our Okinawa data, before I let you off the island.

Next I would like to find the longest trunk road that runs through this prefecture:

SELECT ref, highway, ST_Length(way) as length
    FROM planet_osm_line 
    WHERE highway = 'trunk'
    ORDER BY length DESC;

We have a new function ST_Length() which simply returns the length, given that the geometry is a linestring or multilinestring. The only index that will be used is our "planet_osm_line_highway_index" B-Tree index to perform our Bitmap Index Scan.

ST_Length() does obviously not work with bounding boxes and therefor cannot use the geometrical GiST index. This is yet another function you should use carefully.

When looking at the result set that was returned to us, you will see that some routes show up multiple times. Take route 58, which is the longest and most famous route in Okinawa. It shows up around 769 times. Why?

This is because, especially for a database prepared for mapping, these pieces of geometry are divided over different tiles.

We thus need to accumulate the length of all the linestrings we find that represent pieces of route 58. First, we could try to accomplish this with plain SQL:

WITH road_pieces AS (
    SELECT ST_Length(way) AS length
    FROM planet_osm_line
    WHERE ref = '58'
    )
SELECT sum(length) AS total_length
FROM road_pieces;

This will return:

536468.804010367

Meaning a total length of 536.486 Kilometers. This query will run in about 19.375 ms. Let us add an index to our "ref" column:

CREATE INDEX planet_osm_line_ref_index ON planet_osm_line(ref);

Perform vacuum:

VACUUM ANALYZE planet_osm_line;

This index creation will speed up to query and make it run in little over 3.524 ms. Nice runtime.

You could also perform almost the exact same query, but instead of using an SQL sum() function, you could use ST_Collect(), which creates collections of geometry out of all the separate pieces you feed it. In our case we feed it separate linestrings, which will make this function output a single multilinestring. We would then only have to perform one length calculation.

WITH road_pieces AS (
    SELECT ST_Collect(way) AS geom
    FROM planet_osm_line
    WHERE ref = '58'    
    )
SELECT ST_Length(geom) AS length
    FROM road_pieces;

This query will run even around 1 ms faster then former and it returns the exact same distance of 536.486 Kilometers.

Now that we have this one multilinestring which represents route 58, we could check how close this route comes to our famous Kin building (which we will statically feed):

WITH road_pieces AS (
    SELECT ST_Collect(way) AS geom
    FROM planet_osm_line
    WHERE ref = '58'    
    )
SELECT ST_Distance(geom, ST_GeomFromText('POLYGON((14239931.42 3054117.72,14239990.49 3054224.25,14240230.15 3054091.38,14240171.08 3053984.84,14239931.42 3054117.72))',900913)) AS distance
    FROM road_pieces;

Which would give us:

 7900.58662432767

In other words: Route 58 is, at it closest point, 7.9 Kilometers from our coal storage building. This query now took about 5 ms to complete. A rather nice throughput.

Okay, enough exploring for today.

We took a brief look at indexing our spatial data, and what benefits we could gain from it. And, as you can imagine, a lack of indexes and improper use of the GIS functions, could lead to dramatic slow-downs, certainly on larger data sets.

Shapefiles

Before I will let you go I want to take a brief look at another mechanism of carrying around GIS data: the shapefile. Probably more used then the OSM XML format, but less open. It is almost the GIS standard way of exchanging data between GIS systems.

We can import shapefiles by using a tool called "shp2pgsql" which comes shipped with PostGIS. This tool will attempt to upload ESRI shape data into your PostGIS enables database.

ESRI?

ESRI stands for Environmental Systems Research Institute and is yet another organization that taps into the world of digital cartography.

They have defined a (somewhat open) file format standard that allows the GIS world to save their data in a so called shapefile. These files hold GIS primitives (polygons, linestrings, points, ...) together with a bunch of descriptive information that tells us what each primitive represents.

It was once developed for ESRI's own, proprietary software package (ArcGIS), but was quickly picked up by the rest of the GIS community. Today, almost all serious GIS packages have the ability to read and/or write to such shapefiles.

Shapefile build-up

Let us take a peek at the guts of such a shapefile.

First, contrary to what the name suggest, a shapefile is not a single file. At a minimal level, it is a bundle containing a minimum of three files to be spec compliant:

  • .shp: the first mandatory file has the extension .shp and holds the GIS primitives themselves.
  • .shx: the second important file is an index of the geometry
  • .dbf: the last needed file is a database file with geometry attributes

Getting shapefile data

There are many organizations who offer shapefiles of all areas of the globe, either free or for a small fee. But since we already have data in our database we are familiar with, we could create our own shapefiles.

Exporting with pgsql2shp

Besides "shp2pgsql", which is used to import or load shapefiles, we also got shipped a reverse tool called "pgsql2shp", which can export to or dump shapefiles based on geometry in your database.

So let us, per experiment, create a shapefile containing all secondary roads of Okinawa.

First we need to prepare an empty directory where this tool can dump our data. Since it will create multiple files, it is best to put them in their own spot. Open up a terminal window and go to your favorite directory-making place and create a directory called "okinawa-roads":

$ mkdir okinawa-roads

Next enter that directory.

The "pgsql2shp" tool needs a few parameters to be able to successfully complete. We will be using the following flags:

  • -f, tells the tool which file name to adhere
  • -u, the database user to connect with

After these flags we need to input the database we wish to take a chunk out of and the query which will determine the actual data to be dumped.

The above will result in the following command:

$  pgsql2shp -f secundairy_roads -u postgres gis "select way, ref from planet_osm_line where highway = 'secondary';"

As you can see we construct a query which only gets the road reference and the geometry "way" column from the secondary road types.

After some processing it will have created 4 files, the 3 mandatory ones mentioned above, and a new one called a projection file. This file contains the coordinate system and other projection information in WKT format.

This bundle of 4 files is now our shapefile format which you could easily exchange between GIS aware software packages.

Importing with shp2pgsql

Let us now import these shapefiles back into PostgreSQL and see what happens.

For this we will ignore out "gis" database, and simply create a new database to keep things separated. Connect to a PostgreSQL terminal, create the database and make it PostGIS aware:

CREATE DATABASE gisshape;
\c gisshape
CREATE EXTENSION postgis;

Now go back to your terminal window to do some importing.

The import tool works by dumping the SQL statements to stdin or to a SQL dump file if preferred. If you do not wish to work with such a dump file, you have to pipe the output to the psql command to be able to load in the data.

From the directory where you saved the shapefile dump, run the "shp2pgsql" tool:

$ shp2pgsql -S -s 900913 -I secundairy_roads | psql -U postgres gisshape

Let me go over the flags we used:

  • -S: is used to keep the geometry simple. The tool otherwise will convert all geometry to its MULTI... counterpart
  • -s: is needed to set the correct SRID
  • -I: specifies that we wish the tool to create GiST indexes on the geometry columns

Note that the -S flag will only work if all of your geometry is actual simple and does not contain true MULTI... types of geometry with multiple linestrings, points or polygons in them.

An annoying fact is that you have to tell the loader which SRID your geometry is in. There is a .prj file in our shapefile bundle, but it only contains the WKT projection information, not the SRID. One trick to find the SRID based on the information in the projection file is by using OpenGEO's Prj2EPSG" website, which does quite a good job at looking up the EPSG ID (which most of the time is the SRID). However, it fails to find the SRID of our OSM projection.

Another way of finding our about the SRID is by using the PostGIS spatial_ref_sys table itself:

SELECT srid FROM spatial_ref_sys WHERE srtext = 'PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]';

This will gives us:

900913

Perfect!

If you now connect to your database and query its structure:

\c gisshape
\d

You will see we have a new table called "secondary_roads". This table now holds only the information we dumped into the shapefile, being our road route numbers and their geometry. Neat!

The end

Good.

We are done folks. I hope I have given you enough firepower to be able to commence with your own GIS work, using PostGIS. As I have said in the beginning of this series, the past three chapters form merely an introduction into the capabilities of PostGIS, so as I expect you will do every time: go out and explore!

Try to load in different areas of the world, either with OpenStreetMap or by using shapefiles. Experiment with all the different GIS functions and operators that PostGIS makes available.

And above all, have fun!

And as always...thanks for reading!

by Tim van der Linden at June 25, 2014 10:00 AM

June 20, 2014

Boston GIS (Regina Obe, Leo Hsu)

PostGIS Training sessions at Postgres Open 2014: Chicago

We'll be giving two training sessions on PostGIS at Postgres Open in Chicago September 17th 2014. First will be fairly introductory PostGIS Introduction: Geometry, Geography, and Geocoding and the second will be on more advanced concepts and tools PostGIS on caffeine: Raster, Topology, and pgRouting. Check out the other sessions at Postgres Open Conference Sssions Postgres Open 2014.

by Regina Obe (nospam@example.com) at June 20, 2014 06:50 AM

June 18, 2014

shisaa.jp

Postgis, PostgreSQL's spatial partner - Part 2

Welcome to the secoflynd part of our spatial story. If you have not done so, I advise you to go and read part one first.

The first part of this series gives you some basic knowledge about the GIS world (GIS Objects, WKT, Projections, ...). This knowledge will come in handy in this chapter.

Today we will finally take an actual peek at PostGIS and do some database work:

  • We will see how we can create valid GIS objects and insert them into our database
  • Next let PostGIS retrieve information about these inserted GIS objects
  • Further down the line we will manipulate these object a bit more
  • Then we will leap from geometry into geography
  • Finally we will be doing some real world measurements

Let us get started right away!

Creating the database

Before we can do anything else, we need to make sure that we have the PostGIS extension installed. PostGIS is most of the time packaged as a PostgreSQL contribution package. On a Debian system, it can be installed as follows:

apt-get install postgresql-9.3-postgis-2.1

This will install PostGIS version 2.1 for the PostgreSQL 9.3 database.

Next, fire up your database console and let us first create a new user and database:

CREATE user gis WITH PASSWORD '10gis10';
CREATE DATABASE gis WITH OWNER gis;

Not very original names, I know, but it states its purpose. Next, connect to the gis database and enable the PostGIS extension:

\c gis
CREATE EXTENSION postgis;

Now our database is PostGIS aware, and we are ready to get our hands dirty!

Notice that if you now describe your database:

\d

PostGIS has created a new table and a few new views. This is PostGIS's own bookkeeping and it will store which tables contain geometry or geography columns.

Fun with Polygons

Let us begin this adventure with creating a polygon that has one interior ring, similar to the one we saw in the previous chapter.

Before we can create them, though, we have to create a table that will hold their geometrical data:

CREATE TABLE shapes (
    name VARCHAR
);

Now we have a table named "shapes" with only a column to store its name. But where do we store the geometry?

Because of the new data types that PostGIS introduces (geometry and geography) and to keep its bookkeeping up to date, you can create this column with a PostGIS function named AddGeometryColum():

SELECT AddGeometryColumn('shapes', 'shape', 0, 'POLYGON', 2);

Let us do a breakdown.

First, all the functions that PostGIS makes available to us are divided in groups that define their area of use. AddGeometryColumn() falls in the "Management Functions" group.

It is a function that will create a geometry column in a table of choice and adds a reference to this column to its bookkeeping. It accepts a number of arguments:

  • The table name to where you wish to add the column
  • The actual column name you wish to have
  • The SRID
  • The WKT object you wish to represent
  • The coordinate type you desire (2 means XY)

In the above case we thus wish to add a geometry column to the "shapes" table. The column will be named "shape". The geometry inserted there will get an SRID of 0 and will be of object type POLYGON and have a normal, two dimensional coordinate layout.

SRID?

One thing that you might not yet know from the above function definition is the SRID or Spatial Reference ID and is a very important number when working with spatial data. Remember in the last chapter I kept on yapping about different projections we had and that each projection would yield different results? Well, this is where all this information comes together: the SRID.

Our famous OGC has create a lookup table containing a whopping 3911 entries, each entry with a unique ID, the SRID. This table is called spatial_ref_sys and is, by default, installed into your PostgreSQL database when you enable PostGIS.

But hold on, there is something I neglected to tell you in the previous chapter: the European Petroleum Survey Group or EPSG. The following is something that confuses many people and makes them mix-and-match SRID and EPSG ID's. I will try my best not to add up to that confusion.

EPSG

The EPSG, now called the OGP, is a group of organizations that, among other things, concern themselves over cartography. They are the world's number one authority that defines how spatial coordinates (projected or real world) should be calculated. All the definitions they make get and accompanying ID called the EPSG ID.

The OGC maintains a list to be used inside databases (GIS systems). They give all their entries a unique SRID. These entries refer to defined and official projections, primarily maintained by the EPSG which have their own EPSG ID and unique name. Other projections (not maintained by the EPSG) are also accepted into the OGC SRID list as are your own projections (if you would feel the need).

Let us poke the spatial reference table and see if we can get a more clear picture.

If we would query our table (sorry for the wildcard) and ask for a famous SRID (more on this one later):

SELECT * FROM spatial_ref_sys WHERE srid = 4326;

We would get back one row containing:

  • srid, which is the famous id
  • auth_name, the name of authority organization, in most cases EPSG
  • auth_srid, the EPSG ID the authority organization introduced
  • srtext, tells us how the spatial reference is built using WKT
  • proj4text, commands that drive the proj4 library which is used to make the actual projections

And as you can see, both the "srid" column and the "auth_srid" are identical. This will be the case with many entries.

I should also tell you that this huge list of SRID entries mostly consists of dead or localized projections. Many of the projections listed are not used anymore, but where popular some time in history (they are marked deprecated), or are very localized. In the previous chapter I mentioned that the general UTM system, for example, could be used as a framework for more localized UTM projections. There are hundreds of these local projections that only make sense when used in the area they are intended for.

Simple Features Functions

As I have told you before, the functions that PostGIS makes available are divided into several, defined groups. The functions themselves are too defined, not by PostGIS but by the Simple Features standard maintained by the OGC (as we saw in the previous chapter).

There are a total of 8 major categories available:

  • Management functions: functions which can manipulate the internal bookkeeping of PostGIS
  • Geometry constructors: functions that can create or construct geometry and geography objects
  • Geometry accessors: functions that let us access and ask questions about the GIS objects
  • Geometry editors: functions that let us manipulate GIS objects
  • Geometry outputs: functions that give us various means by which to transform and "export" GIS objects
  • Operators: various SQL operators to query our geography and geometry
  • Spatial relationships and measurements: functions that let us do calculations between different GIS objects
  • Geometry processing: functions to perform basic operations on GIS objects

I have left a few categories out for they are either not part of the Simple Features standard (such as three dimensional manipulations) or beyond the scope. To see a list of all of the functions and their categories, I advise you to visit the PostGIS reference, section 8.

Let us now do some fun manipulations and use some of the functions from these categories, just to get a bit more familiar with how it all works together.

If you inserted the last SQL command which makes the geometry column, you should have gotten back the following result:

public.shapes.shape SRID:0 TYPE:POLYGON DIMS:2

This tells us we created the "shape" column in the "shapes" table and set the SRID to 0. SRID 0 is a convention used to tell a GIS system that you currently do not care about the SRID and simply want to store geometry with an arbitrary X and Y value.

Let us now insert the shape of our square. To insert a polygon into your column, you could use various functions. One of these functions is ST_GeomFromText():

INSERT INTO shapes (name, shape) VALUES (
    'Square with hole',
    ST_GeomFromText('POLYGON ((8 1, 8 8, 1 8, 1 1, 8 1), (6 3, 6 6, 3 6, 3 3, 6 3))', 0)
);

This now inserts our polygon and gives it a name. The ST_GeomFromText() function enables us to enter our polygon object using WKT. This function also accepts a second, optional parameter which is the SRID by which we wish to work. The category of this function is called Geometry Constructors.

You know this polygon has two rings, the exterior and the interior. Let us now ask PostGIS to return only the line that represents the exterior ring:

SELECT ST_ExteriorRing(shape)
    FROM shapes
    WHERE name = 'Square with hole';

And we get back:

0102000020E6100000050000000000000000002040000000000000F03F00000000000020400000000000002040000000000000F03F0000000000002040000000000000F03F000000000000F03F0000000000002040000000000000F03F

Oh my...that is not what we expected. But yet it is correct. This is how PostgreSQL stores geometry/geography. The result is correct, yet unreadable to us humans.

If we wish to get back a readable WKT string, we have to convert it using one of the conversion functions:

SELECT ST_AsText(ST_ExteriorRing(shape))
    FROM shapes
    WHERE name = 'Square with hole';

And we get:

LINESTRING(8 1,8 8,1 8,1 1,8 1)

Aha, that is more like it! This we can read!

We used the ST_ExteriorRing() which falls under the Geometry Accessors category and the ST_AsText() function which resides in the category Geometry Outputs.

Okay, now we wish to know the interior ring:

SELECT ST_AsText(ST_InteriorRingN(shape,1))
    FROM shapes
    WHERE name = 'Square with hole';

The result:

LINESTRING(6 3,6 6,3 6,3 3,6 3)

Notice that the function ST_InteriorRingN() requires you to give the integer of which ring you wish to get, starting from 1. As we have seen before, polygon objects can have multiple interior rings, but only a single exterior one.

Next let us ask all the information about what makes up the shape:

SELECT ST_Summary(shape)
    FROM shapes
    WHERE name = 'Square with hole';

And get back:

Polygon[B] with 2 rings+
  ring 0 has 5 points  +
  ring 1 has 5 points

Ohh, that is pretty cool. We get back a human readable string that explains to us how this particular piece of geometry is build.

Let us now add another polygon:

INSERT INTO shapes (name, shape) VALUES (
    'The intersecting one',
    ST_GeomFromText('POLYGON ((14 1, 15 8, 7 8, 7 1, 14 1))', 0)
);

This is a polygon that will intersect with part of our previous polygon. Let us ask PostGIS if these polygons really intersect each other:

SELECT ST_Intersects(a.shape, b.shape)
    FROM shapes a, shapes b 
    WHERE a.name = 'Square with hole' AND b.name = 'The intersecting one';

If all is well, this will simply return TRUE if they intersect and FALSE if they do not. In this case, it will return TRUE.

The counterpart of our intersect function is ST_Disjoint():

SELECT ST_Disjoint(a.shape, b.shape)
    FROM shapes a, shapes b 
    WHERE a.name = 'Square with hole' AND b.name = 'The intersecting one';

Which will return FALSE in our case.

Let us now add a third polygon which does not intersect our previous two:

INSERT INTO shapes (name, shape) VALUES (
    'The solitary one',
    ST_GeomFromText('POLYGON ((20 20, 20 40, 1 40, 1 20 ,20 20))', 0)
);

This polygon will reside well "above" the other two and does not share any space. Let us now see how far this polygon resides from our first polygon, the one with the hole:

SELECT ST_Distance(a.shape, b.shape)
    FROM shapes a, shapes b 
    WHERE a.name = 'Square with hole' AND b.name = 'The solitary one';

This returns us the number "12", which means they are 12 units apart. And remembering the definition of both shapes, the first shape is 8 units tall and the second shape starts at unit 20. This indeed leaves a gap of 12.

Nice! We have just measured the distance between two objects in a spatial database!

Hmmm, this may mean we are getting closer to knowing the distance to Tokyo...but not yet, we need to play a bit more first.

PostGIS also has the ability to manipulate geometry. Let us, for example, try to move our solitary polygon even further away using the ST_Translate() function under the Geometry Editors category:

SELECT ST_Translate(shape, 0, 10)
    FROM shapes
    WHERE name = 'The solitary one';

The ST_Translate() function will accept the to-be-altered geometry and accepts an X, Y and an optional third dimension.

Running this query will give us a binary representation of a new piece of geometry. The original geometry is not altered. So how can we actually move the geometry that resided in the database?

Simply by using SQL:

UPDATE shapes
    SET shape = ST_Translate(shape, 0, 10)
    WHERE name = 'The solitary one';

Let us now check the new distance:

SELECT ST_Distance(a.shape, b.shape)
    FROM shapes a, shapes b 
    WHERE a.name = 'Square with hole' AND b.name = 'The solitary one';

And get back:

22

Aha! Nice! It has now moved ten units upwards.

Now let us alter the distance once again, but this time we will scale the polygon down:

UPDATE shapes
    SET shape = ST_Scale(shape, 0.5, 0.5)
    WHERE name = 'The solitary one';

If we now check the distance, we will get back:

7

Wow, we have gone from 22 to 7, how did that happen?

Well, it is important to know that the ST_Scale() function currently only supports scaling by multiplying each coordinate. This means that the polygon will not only become smaller or bigger, but will also translate as a result. To know exactly how our new, scaled version of our polygon looks, we can use the ST_Boundary() function which shows us the outer most linestring:

SELECT ST_AsText(ST_Boundary(shape))
    FROM shapes
    WHERE name = 'The solitary one';

And we will get:

LINESTRING(10 15,10 25,0.5 25,0.5 15,10 15)

If we compare that to the same result before scaling (which I handily made ready for you):

LINESTRING(20 30,20 50,1 50,1 30,20 30)

You can see that each value in each coordinate simply was divided by 2. This also clarifies why our polygons are now only 7 units apart (first square stops at 8, the scaled square start at 15).

Okay, okay, I guess we have played enough now. We have seen a small glimpse of the operations you can do on GIS data within PostGIS and seen that PostGIS makes all of this work fairly easy.

I guess we can now take it one step further and start to actually look at some geography!

Fun with the earth

Up until now we have been working with an SRID of 0, which means undefined, inside a geometry column, meaning the data was of type "geometry". Now we want to go out and explore the actual earth, which means we wish to continue in a geographical coordinate system.

This brings us at a crossroad of choices. First, you will need to ask yourself the same question we pondered in chapter one: you wish to work with geometry or geography?

On the one hand we know that geographical measurements are expensive calculations, but most accurate for they are unprojected. On the other hand GIS convention tells us that in any case, we should continue in a geometrical or Cartesian system, simply because...well...it is a convention.

So what do we do?

It all depends on your specific use case.

When working on a "small" scale, say, part of North America, it would make sense to not use geography. Instead, you could (and should) work in a geometrical system using a very accurate projection with SRID 4267 (datum NAD27) or SRID 4269 (datum NAD83) which are both local UTM variants for North America.

Depending on which region you work in, chances are high you have several local projections with their own datum and coordinate system, ready to use. They are very accurate and less expensive to use then direct geography.

For us, however, we will be working on a large scale, for we want to measure a distance that covers much of the globe. You cannot use a local projection or local datum for that.

In such a case you, again, are presented with two options. You could either neglect the convention and simply use geographical data and functions or be nice and adhere to what is agreed upon and work in a Cartesian system.

We will be doing both and we will use the common SRID 4326. This very popular SRID is by heart geographical, for it uses the geographical coordinate system, but can also be used with geometrical data. Confused?

Join the club.

Let me try to clarify.

First, the authority of this SRID is the EPSG and the EPSG ID is identical to the SRID. It uses a popular datum (remember chapter one) called WGS 84 and is referred to as unprojected for it is a geographical representation. This datum is one that is used in GPS systems and is often referred to as a word wide datum.

When you store objects with an SRID of 4326, you are storing them using geographical coordinates aka latitude and longitude. This in contrast to, for example, the former SRID's, like 4267 or 4269, which store their coordinates in UTM values. When you do measurements between two objects carrying this SRID you have two options. You can either do a geographical or a geometrical measurement.

With a geographical measurement there will be no projection and the system will use the WSG 84 datum (the spheroid) to calculate the distance, in three dimensional space. As we have seen before, such a calculation is more expensive and unconventional.

With a geometrical measurement, your geographical coordinates have to be projected on to a flat Cartesian or geometrical plane. This is done automatically when you ask PostGIS to measure distance using one of the more common geometrical functions. When projecting, all GIS systems will use the Plate Carrée projection which means they will use the stored latitude and longitude coordinates directly as an X and Y value.

Let us see this story in action. First we can take a look at the more native geography data. Let us clean our shape table first:

SELECT DropGeometryColumn('shapes', 'shape');

Here we use the DropGeometryColumn() to remove this column from out "shapes" table. Now clear the table:

TRUNCATE shapes;

Next add a new geography column:

ALTER TABLE shapes ADD COLUMN location geography('POINT');

We create a new geography column with the name location in our "shapes" table. We will only be storing Point types.

Notice that the syntax is different and that here we use plain SQL as opposed to the AddGeometryColumn() function from before. Since PostGIS 2 it is possible to create and drop both geometry and geography columns with standard SQL syntax.

If you wish to rewrite our "shape" column addition from the beginning of this chapter, you could write it like this:

ALTER TABLE shapes ADD COLUMN shape geometry('POLYGON',0);

Looks more native and simple, no? Sorry to tell you this so late in the adventure, but now you know the existence of both the functions and the more native SQL syntax. Both will also keep the PostGIS bookkeeping in sync.

Also, for fun, you could do a describe on the table:

\d shapes

And get back:

 Column  |         Type          | Modifiers 
----------+-----------------------+-----------
name     | character varying     | 
location | geography(Point,4326) |

As you can see, the column is of type geography and automatically gets the famous SRID 4326.

Good, let us now try and find an answer to our famous question, How far is Tokyo from my current location. You will be surprised how trivial this will be.

First, as you might suspect, since we are only interested in a point on the earth and not the shape of your location nor Tokyo, we will suffice with a Point object. Next we will need to insert two points into our database, your location and the center of Tokyo, both in geographical coordinates.

Finding Your Location

This means you need to find out your exact latitude and longitude of the place you are at right now.

This could, of course, be done in a myriad of ways: using your cell phone's GPS capabilities, using your dedicated GPS device or using an online map system. I will choose the latter and will be using OpenStreetMap (what else?) to locate my current position.

Open up your favorite web browser and surf to openstreetmap.org. Once there, punch in your address or use the "Where Am I" function. This would give you a point on the map and in the search bar on the left your latitude and longitude coordinate. Take this coordinate and save is as point data into your fresh column:

INSERT INTO shapes VALUES ('My location', ST_GeographyFromText('POINT(127.6791949 26.2124702)'));

The point I am inserting reflects central Naha, the main city of the Okinawa prefecture. Not my current location, but it serves as an illustrative point.

Note that PostGIS expects a longitude as X and latitude as Y. This is many times reversed as what you get back from other sources.

Now you can insert the location of Tokyo, which I conveniently looked up for you:

INSERT INTO shapes VALUES ('Tokyo', ST_GeographyFromText('POINT(139.7530053 35.6823815)'));

Ah, nice! Okay, are you ready to finally, after all the rambling we went through, know the distance? You already know the syntax, punch in the magic:

SELECT ST_Distance(a.location, b.location)
    FROM shapes a, shapes b 
    WHERE a.name = 'My location' AND b.name = 'Tokyo';

In the case you would live in the exact cartographic center of Naha, you will get back:

   st_distance    
------------------
 1557506.28103692

Yeah! This, my lovely folk, is how far you are from Tokyo, at this very moment.

But what is this number you get back?

The result you see here is the distance returned in Meters, meaning, from the point I inserted as "My location", I am 1557506.28 Meters or 1557.50628 Kilometers from Tokyo.

Very neat stuff, would you not say? PostgreSQL just told us how far we are from Tokyo, awesome!

But wait, we are not finished yet. We have now done the most accurate, real geographical distance measurement using expensive geographical calculations.

There is an "in-between" solution before we jump to geometry. PostGIS gives us the ability to replace our spheroid datum with the more classical sphere. The latter has much simpler calculations, but can still return more accurate results them some of the projections.

To redo our calculation from above with a sphere, simply set the spheroid Boolean, a third and optional parameter to the ST_Distance() function, to False:

SELECT ST_Distance(a.location, b.location, False)
    FROM shapes a, shapes b 
    WHERE a.name = 'My location' AND b.name = 'Tokyo';

The result:

   st_distance    
------------------
1557886.68227339

Which is a total distance of 1557.886 Kilometers, a difference of around 300 Meters.

Let us now repeat this story, but use geometry instead. Let us do it the GIS conventional way.

We do not need to recreate our column as a geometry column and insert our data again. We could cheat a little. PostGIS together with PostgreSQL has the unique capability of casting data from one type to another. So without recreating anything, we could simply cast our geography data into geometry on the fly and see what happens.

Note that casting, while very convenient for quick checks, can render an index totally mute. It is therefor important to think ahead and decide if you want to work with geometry or geography, then create the correct column type and use this without casting.

But for our quick and dirty queries, this is fine. Let us continue:

SELECT ST_Distance(a.location::geometry, b.location::geometry)
    FROM shapes a, shapes b 
    WHERE a.name = 'My location' AND b.name = 'Tokyo';

In this query, we cast (::) the geography data inside the "location" columns into geometry. Now we get back:

   st_distance    
------------------
15.3445794209231

Hmm, that is a different result all together. It looks like a much smaller number then before. What is happening?

We just casted our geography to geometry, this means PostGIS will now use a Cartesian system or projection to calculate the distance in a linear way. When using the distance measuring function ST_Distance() on geometry, it will return not meters but the distance expressed in the units the original data was stored in. Since our data is stored with SRID 4326, its units are latitude and longitude. The value you get back is thus degrees.

In the case of the Naha location, this will be 15.344 degrees from Tokyo.

For our human brain this is difficult to imagine, a result in Meters is much more easy to comprehend. So, let us transform this degree value into a metric value.

It is an estimation that one planar degree (in our Cartesian system) equals 111 KM. So the distance now becomes 15.344 degrees times 111: 1703 Kilometers. That is a difference of about 145 Kilometers.

The reason this difference exist is of the projection we are now using. As we have mentioned a few times before, when going from data containing SRID 4326, PostGIS will automatically use the infamous Plate Carrée projection. This projection, as we have seen before, is the least accurate for something like distance measuring.

So let us poke this projection mechanism and try a different, more accurate one, the Lambert, which carries SRID 3587.

To change the projection PostGIS will use, we can use the ST_Transform() function which casts objects to different SRIDs. Note that ST_Transform() only works for geometry objects, so we have to continue to cast our geography location to be able to use them in this function.

SELECT ST_Distance(ST_Transform(a.location::geometry, 3587), ST_Transform(b.location::geometry, 3587))
    FROM shapes a, shapes b 
    WHERE a.name = 'My location' AND b.name = 'Tokyo';

This will gives us:

   st_distance    
------------------
1602392.18109279

Meaning 1602.392 Kilometers, a difference of about 45 Kilometers. That is indeed in between the Plate Carrée and our native geographical measurement.

Another, even more accurate and popular projection is our famous UTM. It can, however, not be used on a world scale. You can only perform measurements within the same UTM zone.

As mentioned in the previous chapter, there are roughly 60 World UTM zones on the earth, but each zone uses their own projection and their own coordinates. This kind of projection is thus not fit for measuring distance on such a large scale.

Let us therefor take this one step further before I leave you to rest. Let us do a measurement with such a UTM projection. We will make a measurement inside of Japan's mainland UTM zone: 54N which has an SRID of 3095.

First we will have to make another point in our database:

INSERT INTO shapes VALUES ('Aomori', ST_GeographyFromText('POINT(140.750616 40.788079)'));

This point represents the city of Aomori in northern Japan, famous for its huge lantern parades.

First let us measure with the native geographical calculations:

SELECT ST_Distance(a.location, b.location)
    FROM shapes a, shapes b 
    WHERE a.name = 'Aomori' AND b.name = 'Tokyo';

This returns:

st_distance    
------------------
573416.203868172

Or 573.416 Kilometers, which is most accurate.

Next, let us throw the good old Plate Carrée projection at it:

SELECT ST_Distance(a.location::geometry, b.location::geometry)
    FROM shapes a, shapes b 
    WHERE a.name = 'Aomori' AND b.name = 'Tokyo';

This will yield

    st_distance    
------------------
5.20224702126502

Which is in degrees again, doing this times 111 Kilometers will yield a total distance of 577.444 Kilometers.

Then let us measure using the correct UTM projection:

SELECT ST_Distance(ST_Transform(a.location::geometry, 3095), ST_Transform(b.location::geometry, 3095))
    FROM shapes a, shapes b 
    WHERE a.name = 'Aomori' AND b.name = 'Tokyo';

This will give us:

  st_distance    
------------------
573228.002047378

Or 573.228 Kilometers and thus only around 200 meters different, in contrast with the Plate Carrée, which was 4 Kilometers different.

You can see that different projections will result in different measurements. It is therefor crucial to know which one to choose. Some are better used on a local scale, like we just did for Japan, others are better on a global scale.

Again, it all comes down to trade-offs and choices.

Okay, yet another big chunk of PostGIS goodness is taken. I suggest a good rest of the mind.

We have seen how we can insert various types of geometry and geography, we saw how to manipulate and question them and we looked at a few real world measurements.

In the next and final chapter, we will be looking at loading some real GIS data from OpenStreetMap into our PostGIS database, take a quick look around my town here in Okinawa and take a deeper look at creating some important indexes.

And as always...thanks for reading!

by Tim van der Linden at June 18, 2014 10:00 AM