Sickels blogg

I would prefer to use postgis, but have to resign to mysql. I want to create a table with spatial information and be able to access it from mapserver. As I am going to store a lot of datapoints, they may be stored just as two columns with x and y values, but as the dataset grows, that starts to be a bit slow.

I started of with a table with more than 20 columns holding point measurements and a lot of metadata related to those measurements. No, it is not practical to split that table up. I tried to add a point column to that table

alter table mobile_rawdata add column coord point

and then I copied the coordinates into that column

update rawdata set coord=point(lon,lat) where not lat is null;

again fine and dandy, but when I tried to put a spatial index on that column, I was told that that was not possible for a column with null values.

So I made a new table

create table mobile_spatial(
OGR_FID serial primary key,
SHAPE geometry not null,
rawdataid integer not null,
addtime timestamp default current_timestamp);

note that the shape column is to be of type geometry, if it is set to another spatial type, ogr will not recognice it. The column names OGR_FID and SHAPE are used as those are default values for ogr. Rawdataid pointing to the id of the old mobile_rawdata table.

Then I tested it with

ogrinfo -so MySQL:radioecology_9,user=user,password=pass mobile_spatial

but ogrinfo still did not understand that this was supposed to be a table with geometry. Reading up a bit more, I found that I needed two more tables in the database: SPATIAL_REF_SYS to store the, well, spatial reference systems and geometry_columns to store information on what kind of geometry the tables had.

I didn’t manage to find a description for those tables but according to gdal.org they will be created when a shape file is imported using ogr2ogr. Trying to import a shape file, I had them created:

CREATE TABLE geometry_columns (
`F_TABLE_CATALOG` varchar(256) DEFAULT NULL,
`F_TABLE_SCHEMA` varchar(256) DEFAULT NULL,
`F_TABLE_NAME` varchar(256) NOT NULL,
`F_GEOMETRY_COLUMN` varchar(256) NOT NULL,
`COORD_DIMENSION` int(11) DEFAULT NULL,
`SRID` int(11) DEFAULT NULL,
`TYPE` varchar(256) NOT NULL
) ENGINE=InnoDB DEFAULT CHARSET=latin1;

and

CREATE TABLE spatial_ref_sys (
`SRID` int(11) NOT NULL,
`AUTH_NAME` varchar(256) DEFAULT NULL,
`AUTH_SRID` int(11) DEFAULT NULL,
`SRTEXT` varchar(2048) DEFAULT NULL
) ENGINE=InnoDB DEFAULT CHARSET=latin1;

the important column in the latter is the SRTEXT, eg for standard WGS84 lat lon, that is

GEOGCS["WGS 84", DATUM["WGS_1984",SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], UNIT["degree", 0.01745329251994328,
AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]]

the SRID field is to be pointed to by the SRID in the geometry_column table.

The geometry_column table has to have the following values:
F_TABLE_NAME : Name of table
F_GEOMETRY_COLUMN: Column in table that holds the geometrical information
COORD_DIMENTION: Dimentionallity of dataset, usually 2
SRID : pointer to the spatial_ref_sys table
TYPE: Type of geometry in the geometry colum, e.g. POINT or POLYGON
the last to columns may be blank.

When learning to use mapserver and mapscript, the first steps have been to display, zoom and pan. Zooming and panning is quite easy, it is just a line like

$map->zoompoint($zoomfactor,$clickpoint,$map->width,$map->height,$mapextent);
where $map is a mapboject, $zoomfactor is a number telling how many times to zoom in or out, positive for zoom in, negative for zoom out, $clickpoint is a pointobject of the location of the click in image coordinates.

Then enter the urge to do a query as well – but not much luck. I reused the same $clickpoint and tried
$map->queryByPoint($clickpoint, MS_SINGLE, 1);
(after doing some adjustments in the map file) but to no luck. Digging into the kmap source code, I understood after a while that the coordinates needed to be in map coordinates and not image coordinates. That is just a simple linear transformation, just with the little snag that the y-coordinates go the opposite direction in the image than in the map. Kmap called the transformation function twice, once for each koordinate, I preferred to make one function for the two-dimentional transformation:

function pix2geo($mapx,$mapy,$minx,$miny,$maxx,$maxy,$sizex,$sizey)

// $mapx and $mapy: image coordinates of point of interest
// $minx, $maxx,$miny, $maxy map extent in map coordinates
// $sizex and $sizey, width and height of image
{
	list($miny,$maxy)=array($maxy,$miny);
	// swaps miny and max y since the coordinates acsends upwards in the map, downwards in the picture
	$geopoint=array();
	foreach( array('x','y') as $dim){
		$p=array();
		foreach(array('map','min','max','size') as $par){
			$name="$par$dim";
			$p[$par]=$$name;
		}
		$geowidth=$p['max']-$p['min'];
		$geopoint[$dim]=$p['map']/$p['size']*$geowidth+$p['min'];
	}
	return(array($geopoint['x'],$geopoint['y']));
}

This can be used like:

list($geoX,$geoY)=pix2geo($_GET['MAP_x']*1,$_GET['MAP_y']*1,$minx,$miny,$maxx,$maxy,$map->width,$map->height);
$geopoint=ms_newpointObj();
$geopoint->setXY( $geoX,$geoY);
$map->queryByPoint($geopoint, MS_SINGLE, 1);

followed by getNumResults for each layer, getResult, getShape and finally getting the relevant values out of the shape.

I have to admit, I simplified things a bit at the beginning. The point needs to be in the projection of the queried data set. As I have the queried dataset in geographical projection and want to have the map in utm, I still have a problem there..

A little phpsnippet for filtering out items from a layer. I have a layer “mobil” which may be on or off. This layer consists of a number of dataset collected at various times. There is one field ‘visible’ that is used to indicate if the dataset is of general interest (eg testing and calibration datasets are set with visible to false) The datasets do also belong to a project, indicated with a projectid in the table. There are also going to be some more items to filter on. Since I don’t know a priori how many items there will be filtered on, I am putting the filter expressions in an array and combining them with ‘and’ before I set the filter on the layer.

if($_GET['mobil']=='ON'){ // no reason to do this if the layer is not visible
	if(!($_GET['invisible']=='ON')){$filters[]="visible = 1";}
	if($_GET['project']*1){$filters[]="projectid = ".$_GET['project']*1;}
	if(count($filters)){ // if there are any items to filter on
		$qLayer = $map->getlayerbyname('mobil');
		$filter="where ".implode(" and ",$filters);
		$qLayer->setfilter($filter);
	}
}

This works for sql-based layers and maybe some others, for some types of layers, one needs to set filter and filteritem.