I’ve been a dormant user of OpenLayers for months (4 months?) now and it was a surprise that the svn trunk had huge differences from what I remember OL (2.4/5?) to be. One of the cool features that and the OpenLayers community contributed was the Google Speherical Mercator hack. Below is a quick step tutorial on how I was able to overlay a custom WMS to Google (set as the baselayer). For this tutorial, I want to overlay a road layer on top of Google.
1. We need to convert our data to Google Projection (Spatial Reference System: 900913). This applies to whatever kind of data (mine is vector stored both in Mapinfo and PostGis) we have. For PostGis, we need to:
INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) VALUES ( 900913,
'spatialreference.org', 900913, '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0
+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs', 'PROJCS
["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6378137,0]],PRIMEM
["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Mercator_2SP"],PARAMETER
["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER
["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],EXTENSION
["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs"]]');
SELECT AddGeometryColumn('public','roads','the_geom_google',900913,'LINESTRING',2);
UPDATE roads SET the_geom_google = Transform(the_geom, 900913);2. MapFile Settings courtesy of SpatialReference: Google Projection
WEB
#Other Web Config Settings goes here...
"wms_srs" "EPSG:900913"
END
PROJECTION
"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
END3. By ensuring that Mapserver has the new 900913 projection, problems such as “msWMSLoadGetMapParams(): WMS server error. Invalid SRS given : SRS must be valid for all requested layers.” or “msProcessProjection(): Projection library error. no options found in ‘init’ file” will be avoided.
cd /ms4w/proj/nad/
gvim epsg
# GCS Voirol 1875 Degree
<104304> +proj=longlat +a=6378249.2 +b=6356514.999904194 no_defs <>
# GCS Voirol Unifie 1960 Degree
<104305> +proj=longlat +ellps=clrk80 no_defs <>
# Google Spherical Mercator
. . .
<900913> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
4. Below is an example WMS Request. Note: “SRS=EPSG 900913″ is added; TRANSPARENT=true not TRANSPARENT=on; Check your BBOX settings for the correct extent.
http://127.0.0.1/cgi-bin/mapserv?
LAYERS=beijing_all
&MAP=%2Fhome%2Fmap%2Fbeijing%2Fnew%2Fbeijing_google.map
&FORMAT=AGG
&TRANSPARENT=true
&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap
&STYLES=
&EXCEPTIONS=application%2Fvnd.ogc.se_inimage
&SRS=EPSG%3A900913
&BBOX=12956687.788758555,4852222.554861524,12956993.536871642,4852528.30297461
&WIDTH=256&HEIGHT=256
5. Requesting the WMS from OpenLayers.
var options = {
projection: "EPSG:900913",
units: "m",
//maxResolution: 156543.0339,
numZoomLevels: 18,
maxExtent: new OpenLayers.Bounds(12823075.86334, 4800551.12375, 13101918.14248, 021301.26141)
};
// avoid pink tiles
OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3;
OpenLayers.Util.onImageLoadErrorColor = "transparent";
map = new OpenLayers.Map('mapdiv',options);
sat_wms = new OpenLayers.Layer.Google(
"Layer",
{type: G_SATELLITE_MAP,'sphericalMercator': true}
);
map.addLayer(sat_wms);
// start custom layer here
var layer_obj = new OpenLayers.Layer.WMS(
"Beijing",
"http://127.0.0.1/cgi-bin/mapserv",
{
layers: 'beijing_all',
map: '/home/map/beijing/new/beijing_google.map',
format: 'AGG',
'transparent': 'true'
}
);
layer_obj.setIsBaseLayer(false);
layer_obj.setVisibility(true);
map.addLayer(layer_obj);
map.addControl(new OpenLayers.Control.MousePosition());
map.addControl(new OpenLayers.Control.LayerSwitcher());
map.addControl(new OpenLayers.Control.Scale());
var center = new OpenLayers.LonLat(12956625.68367, 4852316.90682);
map.setCenter(center, 17);
Here is a couple of Geometric Algorithms used in GIS.
- Convex hull problem: for a set of points, determine the smallest convex set that contains all.
- Line segment intersection: for a set of line segments, determine all intersections.
- Voronoi diagram computation: for a set of points, determine the subdivision of the plane into cells such that inside some cell, one and the same point of the set is closest.
- Delaunay triangulation: for a set of points, determine a planar subdivision by creating edges between the input points in such a way that no two edges intersect, all faces are triangles, no more edges can be added with the given constraints, and no circumcircle of any triangle contains an input point in its interior.
- Minkowski sum: for two simple polygons P and Q, compute the shape that consist exactly of the sum of all points of P and all points of Q, where sum is interpreted as the vector sum.
- Rectangular range search: for a set of points in the plane, design a data structure on those points, such that for every axis-parallel query rectangle, all points in the data structure that lie in the query rectangle can be reported efficiently. Algorithms are needed for the construction of the data structure and for the execution of a query.
Reference:
M Kreveld, Computational Geometry: Its objectives and relation to GIS, Institute of Information and Computing Sciences, Utrecht University
‘R’ is a statistical package. For an overview, please go to www.r-project.org
My intention was to remove the point outliers from a given set of point geometries.
I just recently installed R both on my Windows XP and Debian. Regina’s www.bostongis.com is an excellent tutorial in getting involved with R. I do suggest you head first to PLR Part 1: Up and Running with PL/R (PLR) in PostgreSQL: An almost Idiot’s Guide to get you started.
The install instructions for Windows works flawlessly. I have to hold back to R-2.5 though as I plan to use RPy (Python for R), see details below. To install ‘R’ in Debian, there’s a couple of settings that we need to take care of…
1. Install r-base
sudo apt-get install r-base
2. Install plr on postgres
sudo apt-get install postgresql-8.2-plr
3. Using R in a database
psql -d beijing -U lbs -h 127.0.0.1 < /usr/share/postgresql/8.2/plr.sql
4. Set the R_HOME environment variable
/etc/postgresql/8.2/main/environment
R_HOME='/usr/lib/R'
5. Restart Debian.
RPy, R for Python, is another alternative to use R in Python. I installed it both in Windows and Debian. Note that I reverted to R-2.5 on Windows to be compatible with RPy. For Debian, Im currently using R-2.6.
For the Windows Binary Installation,
1. Read the RPy Main Site
2. Install prerequisites:
- NumPy
- Win32 Extensions Download
3. Afterwards, install the main package, RPy Download
In Debian, its a straight forward…sudo apt-get install python-rpy
I had a problem when using ogr2ogr and converting a postgres table to a road table. My postgres table containa a utf-8 road name which is all in chinese. The mapinfo road table created by ogr2ogr seems to contain the correct geometry and other fields that is in utf-8. However, all my chinese characters is all messed up. So, I have to export the file and open it to mapinfo.
1. In Postgres, to export to a file..
cybersoftbjv1=# set client_encoding = gbk;
SET
cybersoftbjv1=# \o road.txt;
cybersoftbjv1=# select rd_id, cn_name from roads where cn_name <> '';
cybersoftbjv1=# \q
2. Open the file in vim, and do a “%s/ //g”. This would replace all ” ” to “”.
Note: This is reasonable for chinese since chinese dont have spaces. However english prases and sentences differ.
3. Open the file in mapinfo and replace the other columns using Table -> Update.
If anybody has any other way to specify the client encoding in ogr2ogr that would be perfect…
Grabbing the latest ms4w-2.2.6, I was now able to get debugging output by specifying below in my mapfile.
CONFIG MS_ERRORFILE "stderr"
This would log all requests to your Apache error log. I can now play with different debugging modes as described from RFC28. Here is a simple output of my error.log with a DEBUG 2 set..
[error] [client 127.0.0.1] CGI Request 1 on process 4136
[error] [client 127.0.0.1] msDrawMap(): Layer 0 (district), 0.172s
[error] [client 127.0.0.1] msDrawMap(): Layer 1 (water_200k), 0.187s
[error] [client 127.0.0.1] msDrawMap(): Layer 2 (greens_200k), 0.454s
[error] [client 127.0.0.1] msDrawMap(): Layer 6 (roads_150_01), 0.906s
[error] [client 127.0.0.1] msDrawMap(): Layer 9 (roads_270_01), 0.515s
[error] [client 127.0.0.1] msDrawMap(): Layer 10 (roads_180_01), 0.563s
[error] [client 127.0.0.1] msDrawMap(): Layer 11 (roads_280_01), 0.531s
[error] [client 127.0.0.1] msDrawMap(): Layer 12 (roads_400_01), 0.516s
[error] [client 127.0.0.1] msDrawMap(): Layer 13 (roads_140_03), 0.547s
[error] [client 127.0.0.1] msDrawMap(): Layer 14 (roads_140_04), 1.297s
[error] [client 127.0.0.1] msDrawMap(): Layer 15 (roads_141_02), 0.703s
[error] [client 127.0.0.1] msDrawMap(): Layer 16 (roads_150_02), 0.890s
[error] [client 127.0.0.1] msDrawMap(): Layer 17 (roads_boundary_160_170_270), 0.953s
[error] [client 127.0.0.1] msDrawMap(): Layer 18 (roads_160_02), 0.829s
[error] [client 127.0.0.1] msDrawMap(): Layer 19 (roads_170_02), 0.640s
[error] [client 127.0.0.1] msDrawMap(): Layer 20 (roads_270_02), 0.531s
[error] [client 127.0.0.1] msDrawMap(): Layer 22 (roads_180_02), 0.563s
[error] [client 127.0.0.1] msDrawMap(): Layer 23 (roads_280_02), 0.516s
[error] [client 127.0.0.1] msDrawMap(): Layer 24 (roads_400_02), 0.531s
[error] [client 127.0.0.1] msDrawMap(): Layer 26 (subway), 0.109s
[error] [client 127.0.0.1] msDrawMap(): Layer 27 (subwaystops), 0.125s
[error] [client 127.0.0.1] msDrawMap(): Layer 28 (subway_transfer_stops), 0.110s
[error] [client 127.0.0.1] msDrawMap(): Layer 29 (district_boundary), 0.156s
[error] [client 127.0.0.1] msDrawMap(): Layer 33 (400_280_180_170_160_150_labels_01), 4.859s
[error] [client 127.0.0.1] msDrawMap(): Layer 35 (subwaystops_labels), 0.125s
[error] [client 127.0.0.1] msDrawMap(): Layer 36 (district_labels), 0.125s
[error] [client 127.0.0.1] msDrawMap(): Drawing Label Cache, 5.860s
[error] [client 127.0.0.1] msDrawMap() total time: 23.329s