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!

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).

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

Example: Create a 2D linestring column

SELECT AddGeometryColumn ('public', 'segments', 'geom', 2100, 'LINESTRING', 2);

Example: Delete the geometry columns

SELECT DropGeometryColumn('nodes', 'geom');
SELECT DropGeometryColumn('segments', 'geom');

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

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

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

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

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.

UPDATE cities SET geom = ST_SetSRID(ST_MakePoint(lon, lat), 4326);

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.

UPDATE cities SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(lon, lat), 4326), 2100)

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.

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

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.

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

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.

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

There are multiple queries that you can use to solve this problem. You can find an interesting discussion in 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:

SELECT polygons.id FROM polygons WHERE NOT EXISTS
(SELECT 1 FROM points WHERE ST_Intersects(polygons.geom, points.geom))

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.

Example:

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

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.

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

Similar when a CRS transformation is required:

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