一步一步手写GIS开源项目-（1）500行代码实现基础GIS展示功能

1.开篇

       大学毕业工作已经两年了，上学那会就很想研读一份开源GIS的源码，苦于自己知识和理解有限，而市面上也没有什么由浅入深讲解开源gis原理的书籍，大多都是开源项目简介以及项目的简单应用。对于初级程序员想读懂一个成熟的GIS开源项目的困难点主要有三点，1.开发经验和gis原理理解不足。2.一个开源项目是一个循序渐进的过程，如果不是从项目很小的时候跟进，等项目持续更新几年后逻辑就会变得很复杂，小白很难通过Dbug调试来理清楚整个项目的原理。3.GIS属于小行业，国内基本没有国人主导的开源GIS项目，这方面的文章书籍比较欠缺。所以我的想法是通过重写一个GIS项目，还原到项目的初始状态，由浅入深人讲解GIS原理，对于GIS开发从业者不仅要知其然，更要知其所以然。了解GIS底层的开发技术提升自己的竞争力。

     我的文章思路，代码上传到我的github上面，github地址：https://github.com/HuHongYong/ATtuingMap，欢迎大家star 一下，每一篇文章对应的代码生成released版本，方便后期找到文章对应的版本版本代码，如下：

2.开源gis项目的选择

      本次开源项目的选择是SharpMap，选择这个项目的原因主要有两点，1.SharpMap是一套简单易用的小型GIS平台，核心代码1万多行，可以用于开发桌面GIS应用和简单的BS程序。支持多种GIS数据格式，支持空间查询，可输出精美的地图。可了解GIS核心原理。2.开发语言.net，本人工作中使用的核心开发语言，而且现在.net core 3.0开始支持winform、wpf等windows桌面开发技术的跨平台应用的开发，当然这个讲原理就不使用.net core 了，等.net core桌面开发稳定以后可以做一下跨平台。

3.开源GIS最简单实现

     本节500行编写了GIS基本框架，开发环境是vs2017，使用.net framework4.5 ，实现了读取shp点数据，并实现点数据的渲染，以及屏幕坐标向空间坐标转换的基础GIS功能。项目的核心结构如下：

1.Map  只包含一个Map类，是该GIS项目的核心类。

2.Geometries 所有点、线、面等几何类的定义和几何类的方法，这些类都继承了抽象类Geometry，有利于与几何类的扩展，移植，复用。

3.Layers 地图的图层类，该类的核心成员就是下面的4、5、6，构成了一个图层对象的整体。

4.Rendering 提供了用于绘制空间数据的功能，使用c#System.Drawing.Graphics进行渲染绘制。

5.Styles  提供图层的样式，例如点的大小、颜色等。

6.Data 数据的读取接口。例如shapefile文件的读取。

7.Utilities 工具类，一些公用的方法。

4.shp点文件的解析与读取 

       shape文件由ESRI开发，一个ESRI（Environmental Systems Research Institute)的shape文件包括一个主文件，一个索引文件，和一个dBASE表。其中主文件的后缀就是.shp。

      .shp文件包含文件头，数据记录两部分。文件头是固定的100个字节组成，结构如下，数据记录包含着坐标记录信息是文件的数据核心。

       .shp文件头解析，使用的是c# BinaryReader读取字节流，brShapeFile.BaseStream.Seek(36, 0)这个方法是指定位到第36个字节，接着brShapeFile.ReadInt32()是指从第37个字节开始读4个字节。

/// <summary>
/// 读取和解析.shp文件的文件头
/// </summary>
private void ParseHeader (){
     fsShapeFile = new FileStream(_Filename, FileMode.Open, FileAccess.Read);
     brShapeFile = new BinaryReader(fsShapeFile, Encoding.Unicode);
     brShapeFile.BaseStream.Seek(0, 0);
     //读取四个字节，检查文件头
     if (brShapeFile.ReadInt32() != 170328064)
     {
         //文件真实的编码是9994，
         //170328064的16进制为0x0a27，交换字节顺序后就是0x270a，十进制就是9994了
         throw (new ApplicationException("无效的Shapefile文件 (.shp)"));
     }
     //五个没有被使用的int32整数
     brShapeFile.BaseStream.Seek(24, 0);
     //获取文件长度，包括文件头
     _FileSize = 2 * SwapByteOrder(brShapeFile.ReadInt32());
     //读取几何类型
     _ShapeType = (ShapeTypeEnum)brShapeFile.ReadInt32();
     //读取数据的外包矩形
     brShapeFile.BaseStream.Seek(36, 0);
     _Envelope = new BoundingBox(brShapeFile.ReadDouble(), brShapeFile.ReadDouble(), brShapeFile.ReadDouble(),
                     brShapeFile.ReadDouble());

    //通过.shp文件获取数据条数
     // 跳过文件头读取
     brShapeFile.BaseStream.Seek(100, 0);
     // 几何数据记录开始位置
     long offset = 100;

    //遍历数据建立功能包含在数据文件的数量
     while (offset < _FileSize)
     {
         ++_FeatureCount;

        brShapeFile.BaseStream.Seek(offset + 4, 0); //跳过长度
         int data_length = 2 * SwapByteOrder(brShapeFile.ReadInt32());

        if ((offset + data_length) > _FileSize)
         {
             --_FeatureCount;
         }

        offset += data_length; // 添加记录数据长度
         offset += 8; //  +添加每条数据记录头的大小
     }
     _OffsetOfRecord = new int[_FeatureCount];
     //brShapeFile.Close();
     //fsShapeFile.Close();
}

       读取数据记录，由于本次不读取.shx索引文件，我们读取数据生成一个索引数组，方便我们读取数据，通过索引值来读取对应得点数据，代码如下：

/// <summary>
/// 生成矢量文件索引
/// </summary>
private void PopulateIndexes() {
     //记录当前位置的指针
     long old_position = brShapeFile.BaseStream.Position;
     //跳过文件头
     brShapeFile.BaseStream.Seek(100, 0);
     //矢量文件记录开始位置
     long offset = 100;
     for (int x = 0; x < _FeatureCount; ++x)
     {
         _OffsetOfRecord[x] = (int)offset;

        brShapeFile.BaseStream.Seek(offset + 4, 0); //跳过的长度
         int data_length = 2 * SwapByteOrder(brShapeFile.ReadInt32());
         offset += data_length; // 添加记录数据长度
         offset += 8; //   +添加每条数据记录头的大小
     }

    // 返回指针的原始位置
     brShapeFile.BaseStream.Seek(old_position, 0);

}

/// <summary>
/// 从.shp文件中读取并解析几何对象
/// </summary>
/// <param name="oid"></param>
/// <returns></returns>
private Geometry ReadGeometry(int oid) {
     brShapeFile.BaseStream.Seek(_OffsetOfRecord[oid] + 8, 0);
     ShapeTypeEnum type = (ShapeTypeEnum)brShapeFile.ReadInt32(); //Shape type
     if (type== ShapeTypeEnum.Null)
     {
         return null;
     }
     if (type==ShapeTypeEnum.Point)
     {
         return new Point(brShapeFile.ReadDouble(), brShapeFile.ReadDouble());
     }
     else
     {
         throw (new ApplicationException("Shapefile 文件类型 " + _ShapeType.ToString() + " 不支持"));
     }

}

5.屏幕坐标与空间坐标转换

      这一节主要讲一下如何全图展现所有数据。全图显示存在两种状态

1.空间坐标外包矩形宽高比(宽/高)>画布屏幕坐标宽高比(宽/高) 如下图

该状态下展示就以坐标点左右方向占满整个宽度，真实坐标点转为屏幕坐标点的函数。

/// <summary>
/// 空间坐标转屏幕坐标
/// </summary>
/// <param name="p"></param>
/// <param name="map"></param>
/// <returns></returns>
public static PointF WorldtoMap(Point p, Map map)
{
     PointF result = new System.Drawing.Point();

    //在该种情况下，求出的height值为，除去上下空白坐标的高度，如上图所示
     double height = (map.Zoom * map.Size.Height) / map.Size.Width;
     double left = map.Center.X - map.Zoom * 0.5;
     double top = map.Center.Y + height * 0.5 * map.PixelAspectRatio;
     result.X = (float)((p.X - left) / map.PixelWidth);
     result.Y = (float)((top - p.Y) / map.PixelHeight);
     return result;
}

2.空间坐标外包矩形宽高比(宽/高)<画布屏幕坐标宽高比(宽/高) 如下图

该状态下展示就以坐标点左右方向占满整个宽度，真实坐标点转为屏幕坐标点的函数。

/// <summary>
/// 空间坐标转屏幕坐标
/// </summary>
/// <param name="p"></param>
/// <param name="map"></param>
/// <returns></returns>
public static PointF WorldtoMap(Point p, Map map)
{
    PointF result = new System.Drawing.Point();

    //在该种情况下，求出的height值为，整个屏幕高度，如上图所示
     double height = (map.Zoom * map.Size.Height) / map.Size.Width;
     double left = map.Center.X - map.Zoom * 0.5;
     double top = map.Center.Y + height * 0.5 * map.PixelAspectRatio;
     result.X = (float)((p.X - left) / map.PixelWidth);
     result.Y = (float)((top - p.Y) / map.PixelHeight);
     return result;
}

6.绘制点坐标

使用c#System.Drawing.Graphics进行渲染绘制点。

      /// <summary>
       /// 在地图上绘制点
       /// </summary>
       public static void DrawPoint(Graphics g, Point point, Brush b, float size,  Map map) {
           if (point == null)
               return;
           PointF pp = Transform.WorldtoMap(point, map);
           Matrix startingTransform = g.Transform;
           float width = size;
           float height = size;
           g.FillEllipse(b, (int)pp.X - width / 2,
                       (int)pp.Y - height / 2 , width, height);
       }

7总结

      第一节简单的讲了一下.shp数据的读取，以及全图情况下空间坐标与屏幕坐标相互转换。当然只讲了核心功能，具体不明白的可以调试代码进行自己探索。下一节主要讲一下GIS平移缩放问题核心功能。

github项目地址：https://github.com/HuHongYong/ATtuingMap 

作者：ATtuing

出处：http://www.cnblogs.com/ATtuing

本文版权归作者和博客园共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接。