用程序读历史，以数据讲故事。
要说不花钱、不写代码做地理数据可视化分析，最佳选择应该没啥争议，就是QGIS（不仅免费开源、而且全平台覆盖）：
https://www.qgis.org/zh-Hans/site/QGIS首页 ▼
不过这里有一个小问题，QGIS支持打开的数据格式或是GeoJSON、或是Shapefile，像CLIWOC提供的IMMA格式数据虽然包含了经纬度坐标信息，却不能直接被打开使用，因此第一步是要做数据格式的转换：
数据格式转换步骤 ▼
简单来说就是通过程序逐行读取IMMA的文件，并逐字段解析需要的信息，然后再通过GDAL库写入成Shapefile文件。说起来GDAL，是另一个宝藏GIS程序库，提供了包括C、C++、Python、Java、C#、Go、Julia、Lua、JavaScript、TypeScript、Perl、PHP、R、Ruby、Rust等几乎所有主流、非主流程序的调用API，这里的数据转换程序可以挑选任意一门你熟悉的程序语言来写：
https://gdal.org/index.htmlGDAL ▼
对应数据的逐行读取解析，一个比较笨的办法是对照着下面这个表格，挑选自己关心的字段按位置截取：
https://stream-ucm.es/CLIWOC/content/CLIWOC15all.htm字段定义 ▼
但如果把浏览器的控制面板打开：
<table>标签中的字段定义 ▼
所有的字段名，开始、结束位置都是定义在<table>标签里面，因此并不需要在代码里面去写死定义字段的开始、结束位置，只需要去解析这个html页面中的<table>标签，然后根据字段名来获取就好了。如果使用的是Java，可以通过jsoup这个库：
https://jsoup.org/jsoup: Java HTML Parser ▼
解析的时候可以不用关注使用什么字段，把字段的定义全都一股脑加载上：
private Map<String, Field> fields = new HashMap<>();......private void initByHtml(String formatHtml) {Document doc = Jsoup.parse(new File(formatHtml));Element table = doc.select("body > table > tbody").first();Elements rows = table.select("tr");rows.stream().forEach(row -> {try {Field field = new Field(row);fields.put(field.variable, field);} catch (Exception e) {}});}
同样，针对每一行的解析其实也不用太过关注用到什么字段，全部都以key-value的键值对形式存入到一个map中：
public Map<String, String> parseLine(String line) {Map<String, String> values = new HashMap<>();fields.entrySet().stream().forEach(entry -> {Field field = fields.get(entry.getKey());if ("Not in use".equals(field.description)) {return;}try {String value = line.substring(field.start - 1, field.end).trim();values.put(field.variable, value);} catch (Exception e) {}});return values;}
等所有数据解析加载完成，其实就是根据定义好的属性往Shapefile文件里面写：
// 创建图层Layer layer = ds.CreateLayer("layer", srs, ogr.wkbLineString, new Vector<Object>());// 创建属性表// yearFieldDefn year = new FieldDefn("year", ogr.OFTInteger);layer.CreateField(year);// shipIdFieldDefn shipId = new FieldDefn("shipId", ogr.OFTString);shipId.SetWidth(8);layer.CreateField(shipId);// country nameFieldDefn country = new FieldDefn("country", ogr.OFTString);country.SetWidth(16);layer.CreateField(country);
不过等你兴冲冲生成好Shapefile、并用QGIS打开后可能会发现一个问题，就是这些航线这么乱糟糟的甚至还有从陆地上穿过的：
乱糟糟的航线记录 ▼
不用怀疑，一定是程序出了问题，飞机也不带这么直直飞的，何况是200多年前的帆船，太多了可能并不太好分析，咱们可以选择其中的一条船、比如ID为1060的来分析：
ID:1060的航线 ▼
上面这张图是比较能说明问题的，明显中间这根横穿的线条是不应该出现的“幽灵”，原因在于虽然地球是圆的、但GIS软件却并不是，因此当连续2个点穿越地图边缘的时候就出现了不应该有的连线，这里要想正确显示航线就必须在程序里做线条拆分，来看一下坐标点数据：
[[124.58,15.08],[125.15,14.82],[126.2,15.33],[126.71,14.83],[127.78,13.95],[128.23,13.77],[128.66,13.68],[129.4,13.58],[131.03,14.05],[133.56,15.65],[134.66,16.52],[136.08,16.97],[137.16,17.25],[137.45,18.15],[137.58,18.68],[138.41,18.75],[139.45,18.73],[140.6,19.15],[142.3,19.65],[142.55,20.7],[143.71,22.08],[142.98,23.15],[143.5,24.67],[144.45,26],[145.71,27.28],[147.91,29.18],[150.15,30.48],[153,31.82],[153.91,31.37],[154.18,32.38],[156.41,33.38],[158.61,33.98],[160.11,34.43],[160.93,34.85],[164.58,35.37],[167.05,35.88],[169.15,36.28],[171.71,36.4],[172.33,36.8],[171.71,38.07],[172.1,37.7],[172.8,37.62],[173.66,37.4],[174.73,38.13],[176.9,38.6],[179.06,38.7],[180.48,38.62],[183.68,37.87],[187.15,37.1],[190.81,36.25],[194.31,35.23],[197.78,34.85],[200.6,34.4],[203.6,34.77],[206.5,35.23],[210.31,35.7],[212.85,36.1],[214.7,36.18],[217.03,36.55],[218.5,36.42],[219.95,36.65],[221.58,36.95],[222.83,37.22],[224.26,37.38],[226.11,36.85],[227.8,36.08],[229.4,34.9],[232.76,34.9],[234.86,34.98],[235.7,34.68],[235.56,34.73],[238.28,32.53],[238.95,31.45],[239.85,28.93],[241.21,26.73],[242.53,24.85],[244.38,23.85],[245.95,22.82],[247.2,22.35],[249.33,21.75],[250.1,21.55],[251.7,21.13],[253.3,20.3],[254.9,18.43],[255.7,18.2],[255.65,18.08],[256.43,17.83],[257.38,17.33],[258,17],[258.48,16.83]]可以发现CLIWOC采用的经度范围是[0, 360°)，因此线条拆分的时候不仅需要考虑左右边缘、还需要考虑跨越中线的情况：
跨越中线引起的错误 ▼
因此拆分线条的程序逻辑可以这么来写：
List<JSONArray> lines = new ArrayList<>();JSONArray points = new JSONArray();double lastLon = -1;double lastLat = -1;for (Logbook shipLog : shipLogs) {// 获取当前坐标点double lon = shipLog.getLon();double lat = shipLog.getLat();// 第一个点，直接添加if (lastLon < 0) {JSONArray point = new JSONArray();point.add(lon);point.add(lat);points.add(point);lastLon = lon;lastLat = lat;continue;}// 跨越中线无需拆分boolean isCrossMidline = false;if (((lon > 330 && lastLon < 20) // 右到左跨越中线|| (lon < 20 && lastLon > 330)) // 左到右的跨越中线&& Math.abs(lastLat - lat) < 20 // 纬度跨越不能太大) {isCrossMidline = true;}// 拆分线条if (!isCrossMidline&&((lastLon < 180 && lon > 180) // 向右跨越右边界|| (lastLon > 180 && lon < 180) // 向左跨越左边界|| Math.abs(lastLon - lon) > 20 // 单次横向跨越太大|| Math.abs(lastLat - lat) > 20 // 单次纵向跨越太大)) {lines.add(points);// 新航行points = new JSONArray();JSONArray point = new JSONArray();point.add(lon);point.add(lat);points.add(point);lastLon = lon;lastLat = lat;continue;}// 正常添加点至航行JSONArray point = new JSONArray();point.add(lon);point.add(lat);points.add(point);lastLon = lon;lastLat = lat;}lines.add(points);
现在再来生成看一下，是不是航线记录就清爽多了：
最终完成版 ▼
如果仅仅只是做静态地图展示，把前面生成的Shapefile文件用QGIS打开然后设置样式就好，而如果想让航线动起来，我们则还需要往Shapefile里面塞入一个时间：
FieldDefn time = new FieldDefn("time", ogr.OFTDateTime);layer.CreateField(time);
这里先在属性表中定义一个DateTime类型的字段，然后往其中塞入一个ISO 8601格式的日期值：
public String getTime() {StringBuilder sb = new StringBuilder();sb.append(padStart(getYear(), 4, "0")).append("-");sb.append(padStart(getMonth(), 2, "0")).append("-");sb.append(padStart(getDay(), 2, "0")).append("T");sb.append(padStart(getHour(), 2, "0")).append(":");sb.append("00:00Z");return sb.toString();}
处理好后的属性表就像下面这样：
添加了时间字段的属性表 ▼
接下来咱们需要借助于QGIS中的TimeManager这个插件来生成动态效果：
安装TimeManager插件 ▼
下面来看一看动起来后的效果：
大洋上每天的商船位置变化 ▼
跟想象的不一样，别看上面百年航线叠加下大西洋几乎完全被覆盖了，但分散到每一天大洋上的商船并没有几条，如果把时间周期拉长到一年：
大洋上每年的商船航迹 ▼
多的时候也不过是6-7条，所以说数据还是要靠自己动手多维度分析，别人的出来的结果即便是“眼见为实”，或许也并不是看上去的那样。
