通过PROJ4转换地理数据到GoogleMap投影坐标系

Google Map以及VirtualEarth等web gis都采用一种特殊的投影坐标系EPSG:900913，其实这个900913并不是EPSG分配的编号，而是设计Google Map的工程师自己选定的一个编号。该投影坐标系一开始不被EPSG组织承认(EPSG认为这个坐标系的参数设定非常不符合地理科学)，后来因为使用的人越来越多，不得已承认了，但分配了一个别的编号epsg:3785而不是900913。但是大多数程序员不知道，还一直使用900913，呵呵。

关于epsg:3785投影坐标系的详细参数如下：（参考 http://spatialreference.org/ref/epsg/3785/）

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

上面这组参数是用于PROJ4(一个著名的地理投影变换开源库)的，我用以上参数以及从中国国家测绘局下载的中国地图数据和Google Map对照了一下发现x坐标(经度)没有差异，但是y坐标有较大的差异。后来在下述网页http://proj.maptools.org/faq.html的最后一个问题上发现需要一个额外的参数：

 

The coordinate system definition for Virtual Earth Mercator is as follows, which uses a sphere as the earth model for the mercator 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 +no_defs

But, if you do something like:

cs2cs +proj=latlong +datum=WGS84 +to +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 +no_defs

to convert between WGS84 and mercator on the sphere there will be substantial shifts in the Y mercator coordinates. This is because internally cs2cs is having to adjust the lat/long coordinates from being on the sphere to being on the WGS84 datum which has a quite differently shaped ellipsoid.

In this case, and many other cases using spherical projections, the desired approach is to actually treat the lat/long locations on the sphere as if they were on WGS84 without any adjustments when using them for converting to other coordinate systems. The solution is to "trick" PROJ.4 into applying no change to the lat/long values when going to (and through) WGS84. This can be accomplished by asking PROJ to use a null grid shift file for switching from your spherical lat/long coordinates to WGS84.

 

cs2cs +proj=latlong +datum=WGS84 
    +to +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 +no_defs

Note the strategic addition of +nadgrids=@null to the spherical projection definition.

Similar issues apply with many other datasets distributed with projections based on a spherical earth model - such as many NASA datasets, and also (I think) the Google Maps mercator projection.

 

原来缺少+nadgrids=@null参数用来避免sphere lat/long调节就可以了，最终的PROJ4参数如下:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +nagrids=@null +no_defs

这回数据就非常符合Google Map了。