Tuesday, 13 November 2007

Experimenting with MySQL spatial extensions

Two years ago I was experimenting with PostGIS on etch and in these days I'm trying to do the same with MySQL. I'm no longer maintaining a page with my system configuration because It was more annoying than useful (at least for me ;-P ), it should be enough to know that we are on a debian lenny, with MySQL 5.0.45-1.

To test the spatial extensions of MySQL, I'll use again the free GIS data about world borders. The first step is to convert the shape file in a SQL script suitable for MySQL:

  • download shp2mysql from http://kartoweb.itc.nl/RIMapper/ (currently it's at version 0.4);
  • compile it. Assuming that all needed library and header files are installed in the system, it's sufficient to do:
    $ gcc -lshp shp2mysql.c -o shp2mysql
    
  • the usage is analogous to shp2pgsql, for example:
    $ ./shp2mysql -s 4326 ../../geo_data/world_borders/world_borders.shp world_borders db_name > world_borders.mysql.sql
    

Now you you can execute the obtained SQL script in a MySQL database and then query the world_borders table:

> select CNTRY_NAME from world_borders where Contains(ogc_geom, GeomFromText('Point(13 40)')) = 1 ;
+------------+
| CNTRY_NAME |
+------------+
| Italy      | 
+------------+
1 row in set (0.41 sec)

and that's fine.

Without the need of loading GIS data we can try some geometric functions:

> select Glength(GeomFromText('LINESTRING(12 34, 13 35)')) ;
+---------------------------------------------------+
| Glength(GeomFromText('LINESTRING(12 34, 13 35)')) |
+---------------------------------------------------+
|                                   1.4142135623731 | 
+---------------------------------------------------+
1 row in set (0.00 sec)

From MySQL 5.0 reference manual:

In MySQL, the SRID value is just an integer associated with the geometry value. All calculations are done assuming Euclidean (planar) geometry.

Does it mean that you can't do geometric calculation on spheroid with MySQL spatial extensions? I don't think so, maybe it's only a bit more complicated, read this enlighten (for me) thread: MySQL GIS

A final note: both PostgreSQL with PostGIS and MySQL spatial extensions work well under Microsoft Windows. I know it because not every coworker uses GNU/Linux, unfortunately ;-) .

Posted by Nicola Piccinini at 12:18 PM CET in geo/

Sunday, 30 October 2005

Experiments with PostGIS on etch

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 :-( .

Posted by Nicola Piccinini at 6:43 PM CET in geo/

Monday, 3 October 2005

Into the geospatial community

I'm peering into the nice geospatial online community and, some days ago, I ended up installing PostGIS on stakhanov (system configuration).

This quick guide (curiously it is cited even in Web Mapping Illustrated) is a bit out of date but still mostly correct. Only ignore the last lines about sed magic and the relative (broken) link to official PostGIS doc.

Moreover I've found geowanking, a mailing list full of interesting thread. For example "[Geowanking] more google and gis", "[Geowanking] GeoSkating" from which you arrive at an awesome Google Maps mashup and "[Geowanking] googlemaps and WMS".
For coincidence also Paolo Massa has written a post about GeoSkating recently, highlighting the author's Google Maps API experiments.
Application's beauty apart, the interesting fact is the idea to integrate layers on Google Maps from any WMS server: it seems that Google Maps is predisposed for integration with WMS and WFS although this is not well documented presently. This provide a standard way to use a well-done service that, with its extraordinary success, could become a "must have" for commercial web mapping. On Google Maps and GIS web services integration see also:

Posted by Nicola Piccinini at 11:43 AM CEST in geo/