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 http://mappinghacks.com/data/world_borders.zip --23:03:33-- http://mappinghacks.com/data/world_borders.zip => `world_borders.zip' Resolving mappinghacks.com... 216.218.203.219 Connecting to mappinghacks.com|216.218.203.219|: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) - `world_borders.zip' saved [3310260/3310260] pic@stakhanov:~/devel/geodata$ unzip world_borders.zip Archive: world_borders.zip 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 CREATE DATABASE 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 BEGIN [...] COMMIT pic@stakhanov:~/devel/geodata$ psql -d testgis -f /usr/share/postgresql/7.4/contrib/spatial_ref_sys.sql BEGIN [...] COMMIT VACUUM pic@stakhanov:~/devel/geodata$ psql -d testgis -f country.sql BEGIN [...] COMMIT 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)')) ; length2d ----------------- 1.4142135623731 (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 ------------------- 144013.139256991 (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]') ; length2d_spheroid ------------------- 144013.139256991 (1 row) testgis=# /* but changing the latitude does: */ testgis=# select length2d_spheroid(geometryFromText('LINESTRING(12 64, 13 65)'), 'SPHEROID["GRS_1980",6378137,298.257222101]') ; length2d_spheroid ------------------- 121397.522102769 (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)) ; distance ------------------ 144144.751580593 (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)) ; area ------ 1 (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 :-( .