Sunday, 30 October 2005

Experiments with PostGIS on etch

« No more debian-non-US | Main | Languages and Superfluous »

After the upgrading to etch I had to recompile PostGIS (version 1.0.4) for the new version of PostgreSQL:

pic@stakhanov:~$ /usr/lib/postgresql/7.4/bin/postmaster --version
postmaster (PostgreSQL) 7.4.8

I followed the hints in my previous post but had to change something:
$ ./debian/rules config isn't valid any more, one must instead use $ ./debian/rules build and is forced to compile all PostgreSQL's sources.

After this, I tested the new installation with some free GIS data downloaded from Mapping Hacks:

pic@stakhanov:~/devel$ mkdir geodata
pic@stakhanov:~/devel$ cd geodata
pic@stakhanov:~/devel/geodata$ wget
           => `'
Connecting to||:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 3,310,260 (3.2M) [application/zip]

100%[====================================>] 3,310,260    293.94K/s    ETA 00:00

23:03:47 (248.37 KB/s) - `' saved [3310260/3310260]

pic@stakhanov:~/devel/geodata$ unzip
  inflating: world_borders.dbf
  inflating: world_borders.shp
  inflating: world_borders.shx
pic@stakhanov:~/devel/geodata$ /usr/lib/postgresql/7.4/bin/shp2pgsql -s 4326 world_borders.shp country > country.sql
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
pic@stakhanov:~/devel/geodata$ createdb testgis
pic@stakhanov:~/devel/geodata$ createlang plpgsql testgis
pic@stakhanov:~/devel/geodata$ # if the next command works, 
pic@stakhanov:~/devel/geodata$ # PostGIS has been installed correctly
pic@stakhanov:~/devel/geodata$ psql -d testgis -f /usr/share/postgresql/7.4/contrib/lwpostgis.sql
pic@stakhanov:~/devel/geodata$ psql -d testgis -f /usr/share/postgresql/7.4/contrib/spatial_ref_sys.sql
pic@stakhanov:~/devel/geodata$ psql -d testgis -f country.sql
pic@stakhanov:~/devel/geodata$ psql testgis
Welcome to psql 7.4.8, the PostgreSQL interactive terminal.

Type:  \copyright for distribution terms
       \h for help with SQL commands
       \? for help on internal slash commands
       \g or terminate with semicolon to execute query
       \q to quit

testgis=# select cntry_name, area, sum(area(transform(the_geom, 26591))) as calculated_area from country where cntry_name = 'Italy' group by cntry_name, area ;
 cntry_name |  area  | calculated_area
 Italy      | 301230 | 301199476062.85
(1 row)

testgis=# /* note that the "area" column (in square kilometers, included in  
testgis*# shapefile) and the "calculate_area" column (in square meters) are 
testgis*# practically equal */

Other examples, indipendent of the imported data:

testgis=# select length2d(geometryFromText('LINESTRING(12 34, 13 35)')) ;
(1 row)

testgis=# /* that is the usual distance between that two points in a cartesian
testgis*# plan, but: */
testgis=# select length2d_spheroid(geometryFromText('LINESTRING(12 34, 13 35)'), 'SPHEROID["GRS_1980",6378137,298.257222101]') ;  length2d_spheroid
(1 row)

testgis=# /* the result is much bigger because in this case the points' 
testgis*# coordinates are interpreted as longitude, latitude on a big 
testgis*# spheroid.
testgis*# Changing the longitude doesn't affect the result: */
testgis=# select length2d_spheroid(geometryFromText('LINESTRING(42 34, 43 35)'), 'SPHEROID["GRS_1980",6378137,298.257222101]') ;
(1 row)

testgis=# /* but changing the latitude does: */
testgis=# select length2d_spheroid(geometryFromText('LINESTRING(12 64, 13 65)'), 'SPHEROID["GRS_1980",6378137,298.257222101]') ;
(1 row)

testgis=# /* the same length using the distance function and (I hope)
testgis*# an appropriate SRID: */
testgis=# select distance(transform(setSRID(geometryFromText('POINT(12 34)'), 4326), 26591), transform(setSRID(geometryFromText('POINT(13 35)'), 4326), 26591)) ;
(1 row)

Sincerely, I have a lot of doubts regard of the usage of spatial reference identifier (SRID). I chose 4326 because It's suggested in many documents :-P , then, for calculating the Italy's area, I chose 26591 as I saw in spatial_ref_sys table that it is related to that country and has meters as unit of measurement.
Calculating every country's area with this last SRID yields big gaps respect to the values in the shapefile.
Otherwise, using directly 4326 results, for Italy's area, in a small 33.1005615782892:

testgis=# select cntry_name, area, sum(area(the_geom)) as calculated_area from country where cntry_name = 'Italy' group by cntry_name, area ;
 cntry_name |  area  | calculated_area
 Italy      | 301230 | 33.1005615782892
(1 row)

What are these? Degrees?
An experiment shows that, in these conditions, the area is calculated as if the poligons were on a cartesian plan:

testgis=# select area(setSRID(GeometryFromText('POLYGON((12 34, 13 34, 13 35, 12 35, 12 34))'),  4326)) ;
(1 row)

But wasn't 4326 a SRID for longitude, latitude on spheroid? I don't understand.

Anyway, I think that these questions will remain without answer for a long. In fact in this period I am too busy and haven't free time to examine geospatial matters any more :-( .

Posted by Nicola Piccinini at 6:43 PM CET in geo/
2006-Sep-23 CEST: Superfluo
2007-Nov-13 CET: Superfluo