fryktoria.com

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== PostgreSQL - PostGIS ======

===== Introduction =====
This section contains several notes that can be helpful to people who work with GIS systems and/or analyze GIS data.

Although it is highly preferable that your GIS system performs directly on the geometry columns, there are several cases where you need to import numeric data from spreadsheets into latitude-longitude columns and subsequently update the geometry columns manually. These examples will help.  

There are several cases when analyzing GIS data where the coordinates need to appear in a metric system, so transformations need to take place. In the following examples, I use the ubiquitous WGS84 EPSG:4326 as well as the metric GGRS87 EPSG:2100. 

Please note that the examples are not guaranteed to work in all versions/releases of PostgreSQL or PostGIS. Most examples are typical but work at your own risk!  

===== Create and drop coordinate columns =====
When you use database columns where geometry data will be stored, do not create those columns using the normal database commands. **Use the PostGIS commands in order to specify the geometry parameters.**

Examples: Create a 3D point column in EPSG:4326 and a 2D column in GGRS87. You need to specify the database schema ''public'', the table ''nodes'' and ''nodes2'', the geometry columns ''geom'' in both tables, the SRID ''4326'' and ''2100'', the type of layer (points, linestrings, polygons, etc) and the dimensions (2 for x,y and 3 for x,y,z).

<code>
SELECT AddGeometryColumn ('public', 'nodes', 'geom', 4326, 'POINT', 3);
SELECT AddGeometryColumn ('public', 'nodes2', 'geom', 2100, 'POINT', 2);
</code>

Example: Create a 2D linestring column

<code>
SELECT AddGeometryColumn ('public', 'segments', 'geom', 2100, 'LINESTRING', 2);
</code>

Example: Delete the geometry columns

<code>
SELECT DropGeometryColumn('nodes', 'geom');
SELECT DropGeometryColumn('segments', 'geom');
</code>

===== Insert a point with geometry =====

Example: Insert a point in a table with a primary key called ''id'', geometry column ''geom'' and coordinates longitude 23.18 and latitude 38.13

<code>
INSERT INTO nodes (id, geom) VALUES ('point-1', ST_GeomFromText('POINT(23.18 38.13)', 4326));
</code>

Note that the coordinates are in x y format and therefore the first number is ''longitute'' and the second is ''latitude''. Also note that the two numbers are separated with a space.

Example: same as above in 3D, where ''x y z'' are the numeric values to be inserted

<code>
INSERT INTO nodes (id, geom) VALUES ('point-2', ST_GeomFromText('POINT(x y z)', 2100));
</code>

===== Update a WGS84 geom column with data from lat lon columns =====

Suppose that you have in your table two numeric columns ''lat'' and ''lon'' containing WGS84 values and you want to update the ''geom'' column, which was created to store EPSG:4326 geometries.

<code>
UPDATE cities SET geom = ST_SetSRID(ST_MakePoint(lon, lat), 4326);
</code>

===== Update a GGRS87 metric geom column with data from lat lon columns =====

Suppose that you have in your table two numeric columns ''lat'' and ''lon'' containing WGS84 values and you want to update the ''geom'' column, which was created to store EPSG:2100 geometries.

<code>
UPDATE cities SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(lon, lat), 4326), 2100)
</code>

===== Update a geom column with data from metric x,y columns =====

Suppose that you have in your table two numeric columns ''ggrs87_x'' and ''ggrs87_y'' containing GGRS87 values and you want to update the ''geom'' column, which was created to store EPSG:4326 geometries. A transformation is required.

<code>
UPDATE cities SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(ggrs87_x, ggrs87_y), 2100),4326);
</code>

===== Correct a geom after modifying lat, lon =====
Suppose you made a correction to the numeric values of lat and lon of London and you want to update the geom column of this single database row.

<code>
UPDATE cities SET geom = ST_SetSRID(ST_MakePoint(lon, lat), 4326) WHERE city_name = 'London'
</code>

===== Get the nearest neighbours using spatial indexing =====
Suppose that you want to retrieve the 5 nearest neighbours of  point (22.164885,38.924171). The fastest and most effective way is shown below, using the ''<->'' operator, if supported by your versions. 

The data are in EPSG:4326 but it requires the ''ST_SetSRID'' to make sure that both sides of the ''<->'' operator are in the same SRID.
I set both sides to EPSG:2100 because I want to know the distance in meters. Otherwise, I could avoid ST_Transform.

This same query also provides the distance in meters between the provided coordinates and each layer point, as return query column ''dist''. 

<code>
SELECT *, ST_Transform(ST_SetSRID(geom,4326),2100) <-> ST_Transform(ST_SetSRID(ST_MakePoint(22.164885,38.924171),4326),2100) as dist      
FROM cities 
ORDER BY dist
LIMIT 5
</code>

===== Identify the points of a Point layer that are outside polygons of a Polygon layer =====

There are multiple queries that you can use to solve this problem. You can find an interesting discussion in [[https://gis.stackexchange.com/questions/376234/find-polygons-that-has-no-points-in-postgis|Find polygons that has no points in PostGIS]]. For big data sets, I have found that the query that executes in the shortest time is:

<code>
SELECT polygons.id FROM polygons WHERE NOT EXISTS
(SELECT 1 FROM points WHERE ST_Intersects(polygons.geom, points.geom))
</code>

For comparison purposes, in a data set with 45.000 points and 45.000 buildings, this query took approximately 4 minutes, while other queries took more than 14 minutes!

**WARNING:** You need to pay particular attention on how much is exactly "outside". You may encounter cases where the point is considered to be at exactly the perimeter of the polygon, because it has been "snapped" at the side of the polygon, using the snapping features of your GIS software. Yet, there is always a tolerance in snapping and the point might be geometrically within a few tenths of a millimeter inside or outside the polygon. A very small buffer area should be considered around the point to account for this ambiguity. 

===== Import an SHP file manually =====
Example:

<code>
shp2pgsql -I -s <srid> <path> <schema>.<table> -W <encoding>  | psql -d <database>
</code>

===== PHP parameter binding with geometry values =====
If you write PHP code and you want to use the geometry column in a parameter binding, you will need to transform the geometry value into text. 

In the following example, I assume that I have read a value from column ''geom'' with a previous SQL query and now I want to bind the value into parameter '':param'' of a new query.

<code>
$sql = "SELECT ST_AsText(ST_AsEWKT(:param)) AS linecoords";
$st = $dbh->prepare($sql);
$st->bindParam(":param", geom);
$st->execute();
</code>

Similar when a CRS transformation is required:

<code>  
$sql = "SELECT ST_AsText(ST_Transform(ST_AsEWKT(:param), 4326)) AS linecoords";  
$st = $dbh->prepare($sql);
$st->bindParam(":param", geom);
$st->execute();
</code>


~~DISQUS~~
fryktoria.com

User Tools

Site Tools

Page Tools
