openlayers

Using TileCache, OpenLayers, Mapserver for Projection 900913

Apr 8th

Posted by rupert in mapserver

No comments

I had a few problems with TileCache the other week which I am eager to blog about, since I knew for sure that later on, I might encounter the same. I don’t have the exact errors with me right now, so I’m jotting this down from my head…

  1. Classic Resolutions problem. Use extent_type=loose
  2. Can not set image type

UPDATED (JAN 11, 2010): Classic Resolutions problem:

How are resolutions calculated? Assuming we have:

Original:
Lower Left (LL) or minx, miny: 12453557, -5434940
Upper Right (UR) or max, maxy: 16980842, -1180729

maxResolution = (max – minx)/tilesize = (16980842 – 12453557)/512 = 8842.353
where tilesize = 512.

Therefore, we can set/guess for max so that we have maxResolution as a whole number.

Adjusted:
minx, miny: 12453557, -5434940
maxx, maxy: 16980661, -1180729

gives a maxResolution (whole number) of 8842.

Now, you can use 8842 in both the TileCache.cfg and OpenLayers Javascript. More >

, ,

ExtJS and OpenLayers

Feb 19th

Posted by rupert in ExtJS

No comments

It seems that ExtJS and OpenLayers does not blend well together. One of the bug biters that hit me was the way ExtJS was handling arrays. It would be wise for OpenLayers to be agnostic of these frameworks.

for( var i in blocks ) should be transformed to for( var i = 0; i < blocks.length; i++)

I believe OpenLayers Ticket#1362 closely resembles this bug. Thanks to pgiraud for pointing me to the right direction.

,
Travelsite.ph with Google

Reviving an old web map with Google Maps via OpenLayers

Jan 22nd

Posted by rupert in coldfusion

No comments

Travelsite.ph with Google

An old coworker and I worked on a travel portal for the Philippines called travelsite.ph about 4 years ago. We are now given a task of reviving the old web application and even adding mapping functionalities. Back then, the application was using ColdFusion 4.5 and MySQL 3.

Fingers crossed we dropped the app in a ColdFusion 6/7/8 environment with no changes at all. The app still works! Awesome.. how CF really progressed through the years with backward compatibility. The only changes we made was removing the registration/sign up for a quick demo. I just laughed at the oddities and the no brainer features (pertaining to security) that I made when I was starting out.

The database was also intact and have UTM coordinates. We dropped it to a Debian mysql 5 and works flawlessly since its MyISAM. I had the coordinates exported to lon/lat, so I could directly inject it to OpenLayers/Google. After two hours of fiddling around, I got mapping embedded.. hehe.. courtesy of OpenLayers ofcourse.

Here’s a quick reminder to myself…

A. Google WGS 84 Example.

	window.onload = function() {
		var options = {
					projection: "EPSG:4326",
					numZoomLevels: 19,
					maxExtent: new OpenLayers.Bounds(120.8774, 14.3684, 121.1628, 14.7931)
 
				};
 
		// 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}
		);
 
		map.addLayer(sat_wms);
 
		map.addControl(new OpenLayers.Control.MousePosition();
		map.addControl(new OpenLayers.Control.LayerSwitcher());
		map.addControl(new OpenLayers.Control.Scale());
 
		var center = new OpenLayers.LonLat(121.06504,14.65495);
		map.setCenter(center, 16);
	}

B. Google Mercator Projection

		window.onload = function() {
			var options = {
						projection: "EPSG:900913",
						units: "m",
						maxResolution: 156543.0339,
						numZoomLevels: 19,
						maxExtent: new OpenLayers.Bounds(12823075.86334, 4800551.12375, 13101918.14248, 5021301.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, 18);
		}

What’s the difference between both snippets? Obviously projection is one. Since most of my point data is in lon/lat, then the WGS84 example is good if I don’t want to overlay custom precise data. Remember the x shift problem in Google with Openlayers. The Google Mercator example is used when I want to overlay more custom data, particularly polygons/line that needs to fit on Google Layers. For more details, please see my previous blog post.

,
Road Overlay on Google Vector in Forbidden City and Tiananmen, Beijing, China

OpenLayers + Google Spherical Mercator Example

Dec 22nd

Posted by rupert in GIS

No comments

Road Overlay on Google Vector in Forbidden City and Tiananmen, Beijing, China

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"
END

3. 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);
,

Processing Mapinfo Raster JPEG Images using GDAL

Jun 14th

Posted by rupert in GDAL/OGR

No comments

I have a couple of sat images (raw jpegs) from Google that I want to use with Openlayers/Mapserver. The raw jpegs were registered using Mapinfo via GCP (Ground Control Points).

Mapinfo Raster JPEG Images example:

rupert@rupert-winxp /e/home/map/beijing/new/satimages$ ll
-rw-r–r– 1 rupert None 358 Jan 30 03:23 2NE1.TAB
-rw-r–r– 1 rupert None 7.1M Jan 30 02:38 2NE1.jpg
-rw-r–r– 1 rupert None 356 Feb 1 18:56 2NE2a.TAB
-rw-r–r– 1 rupert None 3.8M Feb 1 08:57 2NE2a.jpg

You cannot fully reference 2NE1.TAB as a Mapserver Layer. I tried to use 2NE1.jpg, but the problem its not georeferenced.

rupert@rupert-winxp /e/home/map/beijing/new/satimages
$ gdalinfo 2NE1.jpg
Driver: JPEG/JPEG JFIF
Size is 8650, 6744
Coordinate System is `'
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 6744.0)
Upper Right ( 8650.0,    0.0)
Lower Right ( 8650.0, 6744.0)
Center      ( 4325.0, 3372.0)
Band 1 Block=8650x1 Type=Byte, ColorInterp=Red
Band 2 Block=8650x1 Type=Byte, ColorInterp=Green
Band 3 Block=8650x1 Type=Byte, ColorInterp=Blue

The georeference coordinates of 2NE1.jpg, just like an ESRI World File, is found in 2NE1.TAB…

rupert@rupert-winxp /e/home/map/beijing/new/satimages
$ cat 2NE1.TAB
!table
!version 300
!charset WindowsLatin1
 
Definition Table
  File "2ne1.jpg"
  Type "RASTER"
  (116.38575,39.906105) (349,6619) Label "Pt 1",
  (116.390072,39.93201) (1160,317) Label "Pt 2",
  (116.42786,39.932296) (8210,253) Label "Pt 3",
  (116.4295878,39.90722318) (8522,6358) Label "Pt 4"
  CoordSys Earth Projection 1, 0
  Units "degree"

I found hurting myself in trying to create an ESRI world file from the current MAPINFO Raster TABS. So, I decided to go for GeoTIFF since its native in Mapserver. Using GDAL utilities my only problem is how to put a coordinate system and reference to the raster.

On windows, you could use Frank’s FWTools. For Linux, compile GDAL by source with python would be extremely helpful later on. For installation of GDAL on Linux, we can use Mapserver’s Verbose Installation in Linux Guide.

GDAL – the saviour!.

GDAL utilities is extremely helpful in reprojection, scaling, image mosaics, etc. For now, we will use gdal_translate and gdal_warp. Please RTFM the utilities.

1. Using gdal_translate to specify the gcp’s registered in Mapinfo.

gdal_translate -gcp pixel line easting northing
Add the indicated ground control point to the output dataset. This option may be provided multiple times to provide a set of GCPs.



$ gdal_translate -gcp 349 6619 116.38575 39.906105 -gcp 1160 317 116.390072 39.93201 -gcp 8210 253 116.42786 39.932296 -gcp 8522 6358 116.4295878 39.90722318 -of GTiff 2NE1.jpg 2NE1translated.tif

Input file size is 8650, 6744
0...10...20...30...40...50...60...70...80...90...100 - done.

Note: even if we specify the gcp’s, gdal_translate would not specify the corner coordinates of the tiff.

rupert@rupert-winxp /e/home/map/beijing/new/satimages
$ gdalinfo 2NE1translated.tif
Driver: GTiff/GeoTIFF
Size is 8650, 6744
Coordinate System is `'
GCP Projection =
GCP[  0]: Id=1, Info=
          (349,6619) -> (116.38575,39.906105,0)
GCP[  1]: Id=2, Info=
          (1160,317) -> (116.390072,39.93201,0)
GCP[  2]: Id=3, Info=
          (8210,253) -> (116.42786,39.932296,0)
GCP[  3]: Id=4, Info=
          (8522,6358) -> (116.4295878,39.90722318,0)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 6744.0)
Upper Right ( 8650.0,    0.0)
Lower Right ( 8650.0, 6744.0)
Center      ( 4325.0, 3372.0)
Band 1 Block=8650x1 Type=Byte, ColorInterp=Red
Band 2 Block=8650x1 Type=Byte, ColorInterp=Green
Band 3 Block=8650x1 Type=Byte, ColorInterp=Blue

2. Use gdalwarp to reproject using the gcp and specify the coordinates.

The gdalwarp utility is an image mosaicing, reprojection and warping utility. The program can reproject to any supported projection, and can also apply GCPs stored with the image if the image is “raw” with control information.


$ gdalwarp -s_srs epsg:4326 -t_srs epsg:4326 2NE1translated.tif warped.tif
Creating output file that is 9422P x 5631L.
Processing input file 2NE1translated.tif.
:0...10...20...30...40...5050...60...70...80...90...

Let’s check after gdalwarp using gdalinfo…

$ gdalinfo warped.tif
Driver: GTiff/GeoTIFF
Size is 9422, 5631
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (116.383841930499160,39.933342296341991)
Pixel Size = (0.000004927579869,-0.000004927579869)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( 116.3838419,  39.9333423) (116d23'1.83"E, 39d56'0.03"N)
Lower Left  ( 116.3838419,  39.9055951) (116d23'1.83"E, 39d54'20.14"N)
Upper Right ( 116.4302696,  39.9333423) (116d25'48.97"E, 39d56'0.03"N)
Lower Right ( 116.4302696,  39.9055951) (116d25'48.97"E, 39d54'20.14"N)
Center      ( 116.4070558,  39.9194687) (116d24'25.40"E, 39d55'10.09"N)
Band 1 Block=9422x1 Type=Byte, ColorInterp=Red
Band 2 Block=9422x1 Type=Byte, ColorInterp=Green
Band 3 Block=9422x1 Type=Byte, ColorInterp=Blue

Sweet. Now, all we need to do is display the raster images in Mapserver/OpenLayers.

Specifying a raster image in Mapserver

LAYER
NAME “2NE1″
DATA “satimages/2NE1.tif”
TYPE RASTER
STATUS DEFAULT
END

LAYER
NAME “2NE2″
DATA “satimages/2NE2a.tif”
TYPE RASTER
STATUS DEFAULT
END

Here is the end result…
.

Lessons learned, I tried to specify coordinate extents using gdal_translate -a_ullr ulx uly lrx lry. Specifying the coordinates was subjective by just looking at the cursor location of the registered image in Mapinfo. It is still best to use the GCP’s. Simply put, we need to be accurate in specifying corner coordinates in raster images to project them accurately.

.

, ,
12»