Geoprocessing With Python

0x00. 前言

16 年 12 月研究 GIS 相关资料用于处理 GIS 相关问题,完成基本 GIS 功能。
最新需要进阶相关内容用于更好的处理相关数据。

  • 书籍:
    • Geoprocessing With Python
    • PostGIS In Action 2rd
  • 框架:
    • 前端 Leaflets D3
    • 后端 GeoDjango
  • 其他零零碎碎的资料

特此记录。

本文目录

  • 基本概念
    • Vertor VS Raster
    • Vertor 相关类型与坐标系
    • Raster 相关类型
    • 其他类型
    • GIS 开发的生态圈以及常用技术栈
  • Vertor 分析
  • Raster 分析
  • Vertor 与 Raster

0x01 基本概念

1.1. Vertor VS Raster

  • Vector : 基本单元为 Point : points, lines, and polygons 以及其组合,适用于矢量图,地形边界,路线等。
  • Raster : 基本单元为 Pixel : 2d/3d 包含数值的数组,适用于连续性数据,不仅仅适用于图片。

1.2. Vertor 相关类型与坐标系

1.2.1. 国内常见的几种坐标系

国内由于特殊的国情,国际标准也要向国家标准靠齐。比如各个不同的坐标系上坐标的换算。

我们都知道一个坐标 (x,y) 可以表示为经纬度,甚至放在坐标系上,我们可以这么运算两点 (x1,y1) , (x2,y2) 之间的距离

1
2
# z 表示比例系数
distance = math.sqrt((x1-x2) ** 2 + (y1-y2) ** 2) * z

在近距离的时候的确是可以这么做的比如计算村里小芳和隔壁老王家的距离。当距离过大的时候,比如计算上海 A 区和 B 区的两个写字楼的距离的时候,则有相当大的误差。

那么问题来了:

挖掘技术哪家强

啊不是,是国内有哪些常用坐标标准呢?又是如何计算的呢?

    1. GPS WGS-84 国际标准(原始)
    1. GCJ-02 国内标准(原始数据混淆)
    1. 其他坐标比如 BD-09(原始数据混淆再混淆)

对于小公司而言,我们是没有任何方法来通过 BD-09 以及 GCJ-02 这种坐标系进行运算的:

因为坐标点非线性偏移核心计算方法掌握在 GCJ-02 / BD-09 的公司里面,比如 Google 中国,高德地图,百度地图,腾讯地图。所以,为了研究,则必须要有能够对坐标进行运算的算法, 那这个东西有没有呢?答案是肯定的,因为国外使用的 WGS-84 标准,并且,计算坐标的算法早就开源。

那么,我们的思路就确定下来了。

  1. 各种地图的经纬度坐标比如 BD-09 或 GCJ-02 标转换成 WGS-84 坐标。
  2. 使用开源 GIS 软件进行对 WGS-84 进行运算。

感谢诸多在 GIS 运算上开源的中国先辈,我们轻而易举的获取到了坐标之间相互转化的方法:

https://github.com/wandergis/coordTransform_py

1.2.2. 形状

坐标系,我们可以简单的理解为一个笛卡尔坐标系(虽然这么说很不准确,但已经足够形象了)

于是对于二维的数据,GIS 的分析就可以理解为对于点,线段,多边形自身以及他们之间的关系的分析。

1.3. Raster 相关类型

raster 的 digital elevation model(DEM), 即每一个像素值包含一个 elevation value

GDAL/OGR

0x02. Vertor 分析

0x03. Raster 分析

3.1. 教程

3.2. 教程

3.3. 教程笔记

0x04. Vertor 与 Raster

基础版本

  • Point
  • LineString
  • Polygon

  • MultiPoint

  • MultiLineString
  • MultiPolygon

中级概念

  • Raster / Tile (Bands 是什么鬼)

PostGIS MetaTable

  • spatial_ref_sys
  • geography_columns
  • geometry_columns
  • raster_columns
  • raster_overviews

PostGIS 常用函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
ST_AsText(geom) 用于查看 WKT
ST_GeometryType(geometry) returns the type of the geometry
ST_NDims(geometry) returns the number of dimensions of the geometry
ST_SRID(geometry) returns the spatial reference identifier number of the geometry
ST_X(geometry) returns the X ordinate , 如果作用在 Point 上,则返回经度
ST_Y(geometry) returns the Y ordinate , 如果作用在 Point 上,则返回纬度
ST_Length(geometry) returns the length of the linestring
ST_StartPoint(geometry) returns the first coordinate as a point
ST_EndPoint(geometry) returns the last coordinate as a point
ST_NPoints(geometry) returns the number of coordinates in the linestring
ST_Area(geometry) returns the area of the polygons
ST_NRings(geometry) returns the number of rings (usually 1, more of there are holes)
ST_ExteriorRing(geometry) returns the outer ring as a linestring
ST_InteriorRingN(geometry,n) returns a specified interior ring as a linestring
ST_Perimeter(geometry) returns the length of all the rings
ST_NumGeometries(geometry) returns the number of parts in the collection
ST_GeometryN(geometry,n) returns the specified part
ST_Area(geometry) returns the total area of all polygonal parts
ST_Length(geometry) returns the total length of all linear parts

Snippets

1
2
3
4
5
-- 合并多个区域并返回 multipoly
UPDATE areas as A
SET "Boundary" = ST_Multi(st_union(ARRAY(SELECT geom FROM county_boundary_region WHERE gid in ( 'foo_id','bar_id)')
)))
WHERE A."ID" = 'xxxxxx'

0xEE 参考链接

  1. http://gis.stackexchange.com/questions/6681/what-are-the-pros-and-cons-of-postgis-geography-and-geometry-types
  2. Geo Processing with Python

ChangeLog:

  • 2017-07-11 重修文字