使用1角分高程数据为OpenStreetMap服务器添加海洋等深线

在前文中，我们使用NASA SRTM 数据为陆地添加了3角秒的等高线图层。今天，我们继续为海洋添加等深线。海洋的等深线数据，开放标准的以 etopo1为著名。此数据的分辨率为1角分(arc-min)，即 1/60 度。其dem图层为 WGS-84标准投影，在官方网站可以下载。我选择的是grid配准的binary类型数据。

这个数据其实就是一个16位整形的矩阵，大小为 180*60+1行，360*60+1列，即10801x21601。为了把它变成等深线，导入到postGIS里，需要以下几步。

1 使用QGIS转换投影

WGS-84 的DEM 默认为 EPSG:4326 投影类型，即最通俗的经纬度。在QGIS中打开下载的原始数据，可以看到灰度化的深度图。

我们首先要把它转化为墨卡托投影，使用QGIS控制台命令

gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:3857 -r cubic -wm 1024 -multi -dstnodata -32768 -of GTiff H:/OpenStreetMap/oceanDEM/ETOPO1_Ice_g_gmt4.grd H:/OpenStreetMap/oceanDEM/ETOPO1_Ice_g_gmt4_mkt.tif

或者直接利用图形界面

进行转换。转换后，得到墨卡托投影的高程图。

2 使用QGIS控制台向postGIS计算并导入等高线

这一步走了很多弯路。一开始，我尝试首先生成shp文件，而后利用shp2postgis进行导入。然而，这种方法有2GB大小的文件限制。无奈，经过了一天的研究，发现了超级捷径.

2.1 直接计算、导入等高线

直接利用gdal_contour.exe计算并把结果导入数据库contours，命令：

gdal_contour -a ele -i 10.0 -snodata -32768 D:/GIS/ETOPO1_Ice_g_gmt4.grd -f "postgresql" "PG:host=localhost user=postgres password=XXXXXX dbname=contours"

导入后，生成了一个表 contour

2.2 修改表以适应现有图层模板

为了能够基本照搬上篇文章的SRTM模板，我们把表格稍微修改一下。
1. 修改字段名 id 为 osm_id
2. 删除>=0的元素： delete from contour where ele >=0;
2. 修改geo字段名为way
3. 添加图层开关字段 contour(text), contour_ext
4. 按照不同的高程进行图层分级：

update contour set contour = 'elevation';update contour set contour_ext = 'elevation_minor' where ele % 10 = 0;update contour set contour_ext = 'elevation_medium' where ele % 50 = 0;update contour set contour_ext = 'elevation_major' where ele % 100 = 0;

5.为字段 ele, way 建立索引。

3 添加样式表与图层

与上一章一样，我们为海洋深度添加样式表ocean_depth.mss：

@depth-above: #3A5FCD;@depth-shallow: #7A7A7A;@depth-deep: #778899;@depth-verydeep: #6E6E6E;#depth-minor {    line-width: 0.2;    [elevation >=0] {      line-color: @depth-above;    }    [elevation < 0][elevation > -2000] {      line-color: @depth-shallow;    }    [elevation <=-2000][elevation >-5000] {      line-color: @depth-deep;    }    [elevation <=-5000] {      line-color: @depth-verydeep;    }}#depth-medium {  line-width: 0.4;  [elevation >=0] {    line-color: @depth-above;  }  [elevation < 0][elevation > -2000] {    line-color: @depth-shallow;  }  [elevation <=-2000][elevation >-5000] {    line-color: @depth-deep;  }  [elevation <=-5000] {    line-color: @depth-verydeep;  }}#depth-major {  line-width: 0.6;  [elevation >=0] {    line-color: @depth-above;  }  [elevation < 0][elevation > -2000] {    line-color: @depth-shallow;  }  [elevation <=-2000][elevation >-5000] {    line-color: @depth-deep;  }  [elevation <=-5000] {    line-color: @depth-verydeep;  }}#depth-medium-text {  text-name: "[ele]";  text-face-name: @book-fonts;  text-placement: line;  text-spacing: 500;  text-size: 7;  text-fill: #1a1a1a;  [elevation >=0] {    text-fill: @depth-above;  }  [elevation < 0][elevation > -2000] {    text-fill: @depth-shallow;  }  [elevation <=-2000][elevation >-5000] {    text-fill: @depth-deep;  }  [elevation <=-5000] {    text-fill: @depth-verydeep;  }}#depth-major-text {  text-name: "[ele]";  text-face-name: @book-fonts;  text-placement: line;  text-spacing: 1000;  text-size: 8;  text-fill: #1c1c1c;  [elevation >=0] {    text-fill: @depth-above;  }  [elevation < 0][elevation > -2000] {    text-fill: @depth-shallow;  }  [elevation <=-2000][elevation >-5000] {    text-fill: @depth-deep;  }  [elevation <=-5000] {    text-fill: @depth-verydeep;  }}

而后，在project.mml中引入图层：

...  osm2pgsql: &osm2pgsql    type: "postgis"    dbname: "gis"    key_field: ""    geometry_field: "way"    extent: "-20037508,-20037508,20037508,20037508"  contour: &contour    type: "postgis"    dbname: "contours"    key_field: ""    geometry_field: "way"    extent: "-20037508,-20037508,20037508,20037508"Stylesheet:...  - "contour.mss"  - "ocean_depth.mss"...Layer:  - id: "depth-minor"    name: "depth-minor"    class: ""    geometry: "linestring"    <<: *extents    Datasource:      <<: *contour      table: |-        (SELECT            way,            CAST (ele AS INTEGER) AS elevation,            ele          FROM contour          WHERE contour_ext='elevation_minor' AND ele <= -30          ) AS contours_minor    properties:      minzoom: 14    advanced: {}  - id: "depth-medium"    name: "depth-medium"    class: ""    geometry: "linestring"    <<: *extents    Datasource:      <<: *contour      table: |-        (SELECT            way,            CAST (ele AS INTEGER) AS elevation,            ele          FROM contour          WHERE contour_ext='elevation_medium' AND ele <= -30          ) AS contours_medium    properties:      minzoom: 12    advanced: {}  - id: "depth-major"    name: "depth-major"    class: ""    geometry: "linestring"    <<: *extents    Datasource:      <<: *contour      table: |-        (SELECT            way,            CAST (ele AS INTEGER) AS elevation,            ele          FROM contour          WHERE contour_ext='elevation_major' AND ele <= -30          ) AS contours_major    properties:      minzoom: 8    advanced: {}  - id: "depth-medium-text"    name: "depth-medium-text"    class: ""    geometry: "linestring"    <<: *extents    Datasource:      <<: *contour      table: |-        (SELECT            way,            CAST (ele AS INTEGER) AS elevation,            ele          FROM contour          WHERE contour_ext='elevation_medium'  AND ele <= -30          ) AS contours_medium_text    properties:      minzoom: 13    advanced: {}  - id: "depth-major-text"    name: "depth-major-text"    class: ""    geometry: "linestring"    <<: *extents    Datasource:      <<: *contour      table: |-        (SELECT            way,            CAST (ele AS INTEGER) AS elevation,            ele          FROM contour          WHERE contour_ext='elevation_major' AND ele <= -30          ) AS contours_major_text    properties:      minzoom: 8    advanced: {}  - id: "world" ...

别忘了编译mml为mapnik:

carto ./project.mml > mapnik.xml

4 测试效果

重新渲染后，效果还不错！虽然1角分的分辨率比较低，但是在海洋上还是够用啦！

相关数据与虚拟机镜像在这里下载。
网站测试：
http://www.goldenhawking.org:8088/