基于Python的地理空间分析（三）：GIS数据类型介绍

发表于 分类于 阅读次数： Disqus：

前面已经介绍了GIS数据库和Python的相关库，本文介绍GIS中的主要数据类型，以及如何用Python读写GIS数据。

GIS数据类型

两种数据类型

GIS以数字数据的形式描述现实世界的客观对象，而现实对象可以不严谨地抽象为连续的领域和离散的对象两大类，打个比方，全球的温度、降水是空间里的连续域，而全球的200多个国家和地区则是离散的个体。相应地，针对这两类抽象概念，地理学中也有两种数据模型：

  1. 场（field based）模型，也称为域模型，是把地理空间中的对象作为连续的变量或体来看待，根据应用的不同，场可以表现为二维或者三维，一个二维场就是在二维空间中任意一个位置上，都有一个表征的属性值，即A=f(x,y)。所以场或域，本质上指的是属性在在空间上的连续变化，在实际应用中，真正的“连续”是难以获得和分析的，所以往往通过采样的方法来表征场的变化。
  2. 对象（object based）模型，也称为要素模型，是将研究的地理空间看作一个空域，而地理现象和空间实体作为独立的对象分布其中。按照空间特征，可以分为点、线、面、体四种基本对象，也可以进一步组合构成复杂独享。因此，对象模型适合对具有明确边界的对象进行建模，可以按照空间、时间、非空间属性以及其他要素在时空、语义上的关系来描述。

进一步地，在GIS中对于两种数据模型分别有两大类实现形式：

  1. 栅格数据，适用于场模型抽象。前面提到，域模型需要经过离散化采样才能为我所用，而“栅格”就是采样，栅格单元的形状通常是正方形或者矩形，一个栅格单元如同图像中的一个像素点，组合起来就成为栅格图像，在栅格图像中，可以按照不同的数据进行分层，每层描述一种属性。
  2. 矢量数据，适用于对象模型抽象。矢量数据使用点、线、面等实体及其集合表示地理对象，可以明确描述图形要素的拓扑关系，具有缩放不变性、便于几何计算和占用存储空间小等优点。

简而言之，矢量数据是图形，栅格数据是图像，这就是GIS中的两需要类数据。

所有GIS数据都由一个或另一个组成，但也可以组合使用矢量和栅格。在决定使用哪种数据类型时，请考虑数据所代表的地理信息的规模和类型，这又决定了要使用的Python库。当然，具体选择某个Python库也可能取决于个人偏好，并且可能有多种方法来执行相同的任务。

栅格数据类型

在GIS中，比较流行的栅格数据格式包括(但不限于）以下几种：

  • ECW(Enhanced Compressed Wavelet，增强压缩小波）：ECW是一种压缩图像格式，通常用于航拍和卫星图像。ECW格式以其能够实现高压缩比的同时保持图像对比度质量而闻名。
  • Esri网格（Esri gird）:用于将属性数据添加到栅格文件的文件格式，可以用作整数和浮点网格。
  • GeoTIFF(Geographic Tagged Image File Format，地理标记图像文件格式）：用于GIS和卫星遥感应用的行业图像标准，几乎所有的GIS和图像处理软件包都具有GeoTIFF兼容性。
  • JPEG 2000：一种开源压缩栅格格式（文件扩展名JP2），允许有损压缩和无损压缩，可以实现20:1的压缩比（与MrSID格式相近）。
  • MrSID（Multi-Resolution Seamless Image Database，多分辨率无缝图像数据库）:一种允许无损压缩和有损压缩的压缩小波格式， MrSID图像具有SID的扩展名，并附有文件扩展名为SDW的世界文件。此外，LizardTech专有的MrSID格式通常用于需要压缩的正射影像。

矢量数据类型

Shapefile

Shapefile可能是当今地理矢量数据最常用的数据格式。Shapefile由Esri基于Esri和其他GIS软件产品之间的数据互操作性的开放规范开发。虽然为了替取代hapefile而诞生了许多其他文件格式，但它仍然是一种广泛使用的文件格式。目前，Python中存在许多用于读取和编写shapefile的第三方编程模块。虽然名称shapefile可能暗示只有一个文件与之关联，但是单个shapefile实际上至少需要三个文件需要存储在同一目录中才能正常工作：

  1. 具有要素几何本身的**.shp文件;
  2. 具有要素几何的位置索引的.shx文件；
  3. 允许快速向前和向后搜索具有每个形状的柱状属性的.dbf文件

Shapefile具有自己的结构。主文件（.shp）包含几何数据，由一个固定长度的标头组成，后跟一个或多个可变长度记录。

GeoJSON

GeoJSON是一种基于JSON的文件格式，很快就在GIS领域流行起来。GeoJSON使用JavaScript Object Notation（JSON）开放数据标准将地理要素存储为键值对。正如JSON文件一样，GeoJSON文件具有很好的可读性，可以使用简单的文本编辑器创建，因此在现在的空间数据库、开放数据平台和商业GIS软件中广泛应用。GeoJSON可以用于保存点、线、多边形等各种类型的地理空间矢量数据。此外，GeoJSON不一定使用.geojson作为文件扩展名，还可以使用.json。

KML

KML（Keyhole Markup Language）中的K，指的是开发该格式的公司。KML用于存储地理数据，可以使用Google Earth、ArcGIS、Photoshop和AutoCAD等大量应用程序进行数据的可视化。KML基于XML开发，使用具有嵌套元素和属性的基于标记的结构。 KML文件通常分布在KMZ文件中，这些文件是带有.kmz扩展名的压缩KML文件。 KML的坐标参考系使用由1984年世界大地测量系统（World Geodetic System of 1984, WGS84）定义的经度，纬度和高度坐标。

GeoPackage

GPKG（Open Geospatial Consortium GeoPackage, 开放地理空间联盟GeoPackage）是一种开放数据格式，用于支持矢量和栅格数据的地理信息系统。 该格式由OGC定义并于2014年发布，之后得到了许多政府，商业和开源组织的广泛支持。 GeoPackage数据格式是在考虑移动用户的情况下开发的——它被设计为尽可能高效，所有信息都包含在一个文件中。 这使得在云存储和USB驱动器上快速共享它们变得容易，并且它被用于断开连接的移动应用程序。 GeoPackage文件构建为扩展的SQLite 3数据库文件（* .gpkg），它结合了数据和元数据表。

向量数据读写

使用GeoPandas库

准备工作

最最重要的准备工作就是首先要安装好GeoPandas，为什么要单独拿出来说？如果你是Windows系统，然后没有用虚拟环境，还用pip安装的话，那可真的是太麻烦了！！ 相应地，解决方法有三个，我建议：

用conda安装！

用conda安装！！

用conda安装！！！

除了安装依赖的包以外，如果要下载一些数据的话，推荐Natural Earth，可以找到不同格式、不同分辨率的基本数据。

读写方法

1
2
3
4
5
6
7
8
9
10
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pylab import rcParams

rcParams['figure.figsize'] = 15, 12

%matplotlib inline

1
2
3
df = gpd.read_file(r'data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp', encoding='utf-8')
df.head()
df.plot(figsize=[15,12])

output_20_1

查看GeoPandas数据中的属性，可以用geom_type:

1
df.geom_type.head()
0    MultiPolygon
1    MultiPolygon
2    MultiPolygon
3         Polygon
4    MultiPolygon
dtype: object

查看GeoPandas数据所使用的坐标系，可以用crs(coordinate reference system, CRS)，CRS提供了有关空间数据集的基本信息。：

1
df.crs
{'init': 'epsg:4326'}

上面的结果 epsg:4326 是由国际石油和天然气生产商协会（International Association of Oil and Gas Producers, IOGP）定义的编码，更多信息可以在EPSG网站（www.spatialreference.org ）上查询。EPSG 4326更为人所知的名字叫“WGS 1984”，也就是标准的地球坐标系。

如果想要改变坐标系，可以用to_crs()命令，例如我们改成墨卡托投影坐标系：

1
2
merc = df.to_crs({'init': 'epsg:3395'})
merc.plot(figsize=[15,12])

output_27_1

如果想将读取的shapefile数据转换为json数据，可以采用to_json()命令，to_json()命令只是将数据转换为json格式，如果要保存到文件，则可以用to_file()命令：

1
df.to_json()[0:500]
'{"type": "FeatureCollection", "features": [{"id": "0", "type": "Feature", "properties": {"featurecla": "Admin-0 country", "scalerank": 5, "LABELRANK": 2, "SOVEREIGNT": "Indonesia", "SOV_A3": "IDN", "ADM0_DIF": 0, "LEVEL": 2, "TYPE": "Sovereign country", "ADMIN": "Indonesia", "ADM0_A3": "IDN", "GEOU_DIF": 0, "GEOUNIT": "Indonesia", "GU_A3": "IDN", "SU_DIF": 0, "SUBUNIT": "Indonesia", "SU_A3": "IDN", "BRK_DIFF": 0, "NAME": "Indonesia", "NAME_LONG": "Indonesia", "BRK_A3": "IDN", "BRK_NAME": "Indone'

1
2
3
4
5
6
7
8
9
from shapely import geometry
upcast_dispatch = {geometry.Point: geometry.MultiPoint,
                   geometry.LineString: geometry.MultiLineString,
                   geometry.Polygon: geometry.MultiPolygon}
def maybe_cast_to_multigeometry(geom):
    caster = upcast_dispatch.get(type(geom), lambda x: x[0])
    return caster([geom])
df = df.geometry.apply(maybe_cast_to_multigeometry)
df.to_file(driver='GeoJSON', filename=r'data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.json')
C:\ProgramData\Anaconda3\lib\site-packages\geopandas\io\file.py:108: FionaDeprecationWarning: Use fiona.Env() instead.
  with fiona.drivers():

注：需要增加自定义的mayber_cast_to_multigeometry()函数，将polygon转换为multipolygon。

GeoPandas依赖Fiona库实现文件转换，如果想要查看GeoPandas文件支持的数据转换，也就是Fiona支持的数据格式，采用如下命令：

1
2
import fiona
fiona.supported_drivers
{'AeronavFAA': 'r',
 'ARCGEN': 'r',
 'BNA': 'raw',
 'DXF': 'raw',
 'CSV': 'raw',
 'OpenFileGDB': 'r',
 'ESRI Shapefile': 'raw',
 'GeoJSON': 'rw',
 'GPKG': 'rw',
 'GML': 'raw',
 'GPX': 'raw',
 'GPSTrackMaker': 'raw',
 'Idrisi': 'r',
 'MapInfo File': 'raw',
 'DGN': 'raw',
 'PCIDSK': 'r',
 'S57': 'r',
 'SEGY': 'r',
 'SUA': 'r'}

OGR库读写矢量数据

命令行做法

OGR中使用ogrinfo命令描述矢量数据，但是ogrinfo是命令行而并非python的命令，所以在Jupter Notebook中，在使用前增加感叹号（！）作为前缀。

ogrinfo --formats列出ogrinfo可以访问的可用格式。 结果还告诉我们GDAL / OGR是否只能读取/打开格式，或者它是否也能以该格式写入新图层。 从输出中可以看出，有许多支持OGR的文件格式。 查看列表中的Esri shapefile，添加（rw + v）意味着OGR支持Esri shapefile的读取，写入，更新（意味着创建）和虚拟格式。

1
!ogrinfo --formats
Supported Formats:
  netCDF -raster,vector- (rw+s): Network Common Data Format
  MSSQLSpatial -vector- (rw+): Microsoft SQL Server Spatial Database (BCP)
  PCIDSK -raster,vector- (rw+v): PCIDSK Database File
  JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library
  PDF -raster,vector- (rw+vs): Geospatial PDF
  MBTiles -raster,vector- (rw+v): MBTiles
  EEDA -vector- (ro): Earth Engine Data API
  DB2ODBC -raster,vector- (rw+): IBM DB2 Spatial Database
  ESRI Shapefile -vector- (rw+v): ESRI Shapefile
  MapInfo File -vector- (rw+v): MapInfo File
  UK .NTF -vector- (rov): UK .NTF
  OGR_SDTS -vector- (rov): SDTS
  S57 -vector- (rw+v): IHO S-57 (ENC)
  DGN -vector- (rw+v): Microstation DGN
  OGR_VRT -vector- (rov): VRT - Virtual Datasource
  REC -vector- (ro): EPIInfo .REC
  Memory -vector- (rw+): Memory
  BNA -vector- (rw+v): Atlas BNA
  CSV -vector- (rw+v): Comma Separated Value (.csv)
  NAS -vector- (rov): NAS - ALKIS
  GML -vector- (rw+v): Geography Markup Language (GML)
  GPX -vector- (rw+v): GPX
  KML -vector- (rw+v): Keyhole Markup Language (KML)
  GeoJSON -vector- (rw+v): GeoJSON
  GeoJSONSeq -vector- (rw+v): GeoJSON Sequence
  ESRIJSON -vector- (rov): ESRIJSON
  TopoJSON -vector- (rov): TopoJSON
  OGR_GMT -vector- (rw+v): GMT ASCII Vectors (.gmt)
  GPKG -raster,vector- (rw+vs): GeoPackage
  SQLite -vector- (rw+v): SQLite / Spatialite
  ODBC -vector- (rw+): ODBC
  WAsP -vector- (rw+v): WAsP .map format
  PGeo -vector- (ro): ESRI Personal GeoDatabase
  PostgreSQL -vector- (rw+): PostgreSQL/PostGIS
  OpenFileGDB -vector- (rov): ESRI FileGDB
  XPlane -vector- (rov): X-Plane/Flightgear aeronautical data
  DXF -vector- (rw+v): AutoCAD DXF
  CAD -raster,vector- (rovs): AutoCAD Driver
  Geoconcept -vector- (rw+v): Geoconcept
  GeoRSS -vector- (rw+v): GeoRSS
  GPSTrackMaker -vector- (rw+v): GPSTrackMaker
  VFK -vector- (ro): Czech Cadastral Exchange Data Format
  PGDUMP -vector- (w+v): PostgreSQL SQL dump
  OSM -vector- (rov): OpenStreetMap XML and PBF
  GPSBabel -vector- (rw+): GPSBabel
  SUA -vector- (rov): Tim Newport-Peace's Special Use Airspace Format
  OpenAir -vector- (rov): OpenAir
  OGR_PDS -vector- (rov): Planetary Data Systems TABLE
  WFS -vector- (rov): OGC WFS (Web Feature Service)
  WFS3 -vector- (ro): OGC WFS 3 client (Web Feature Service)
  HTF -vector- (rov): Hydrographic Transfer Vector
  AeronavFAA -vector- (rov): Aeronav FAA
  Geomedia -vector- (ro): Geomedia .mdb
  EDIGEO -vector- (rov): French EDIGEO exchange format
  GFT -vector- (rw+): Google Fusion Tables
  SVG -vector- (rov): Scalable Vector Graphics
  CouchDB -vector- (rw+): CouchDB / GeoCouch
  Cloudant -vector- (rw+): Cloudant / CouchDB
  Idrisi -vector- (rov): Idrisi Vector (.vct)
  ARCGEN -vector- (rov): Arc/Info Generate
  SEGUKOOA -vector- (rov): SEG-P1 / UKOOA P1/90
  SEGY -vector- (rov): SEG-Y
  XLS -vector- (ro): MS Excel format
  ODS -vector- (rw+v): Open Document/ LibreOffice / OpenOffice Spreadsheet
  XLSX -vector- (rw+v): MS Office Open XML spreadsheet
  ElasticSearch -vector- (rw+): Elastic Search
  Walk -vector- (ro): Walk
  Carto -vector- (rw+): Carto
  AmigoCloud -vector- (rw+): AmigoCloud
  SXF -vector- (rov): Storage and eXchange Format
  Selafin -vector- (rw+v): Selafin
  JML -vector- (rw+v): OpenJUMP JML
  PLSCENES -raster,vector- (ro): Planet Labs Scenes API
  CSW -vector- (ro): OGC CSW (Catalog  Service for the Web)
  VDV -vector- (rw+v): VDV-451/VDV-452/INTREST Data Format
  GMLAS -vector- (rwv): Geography Markup Language (GML) driven by application schemas
  MVT -vector- (rw+v): Mapbox Vector Tiles
  TIGER -vector- (rw+v): U.S. Census TIGER/Line
  AVCBin -vector- (rov): Arc/Info Binary Coverage
  AVCE00 -vector- (rov): Arc/Info E00 (ASCII) Coverage
  NGW -raster,vector- (rw+s): NextGIS Web
  HTTP -raster,vector- (ro): HTTP Fetching Wrapper

1
!ogrinfo -so "data/ne_10m_admin_0_countries" ne_10m_admin_0_countries
INFO: Open of `data/ne_10m_admin_0_countries'
      using driver `ESRI Shapefile' successful.

Layer name: ne_10m_admin_0_countries
Metadata:
  DBF_DATE_LAST_UPDATE=2018-05-21
Geometry: Polygon
Feature Count: 255
Extent: (-180.000000, -90.000000) - (180.000000, 83.634101)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
featurecla: String (15.0)
scalerank: Integer (1.0)
LABELRANK: Integer (2.0)
SOVEREIGNT: String (32.0)
SOV_A3: String (3.0)
ADM0_DIF: Integer (1.0)
LEVEL: Integer (1.0)
TYPE: String (17.0)
ADMIN: String (36.0)
ADM0_A3: String (3.0)
GEOU_DIF: Integer (1.0)
GEOUNIT: String (36.0)
GU_A3: String (3.0)
SU_DIF: Integer (1.0)
SUBUNIT: String (36.0)
SU_A3: String (3.0)
BRK_DIFF: Integer (1.0)
NAME: String (25.0)
NAME_LONG: String (36.0)
BRK_A3: String (3.0)
BRK_NAME: String (32.0)
BRK_GROUP: String (17.0)
ABBREV: String (13.0)
POSTAL: String (4.0)
FORMAL_EN: String (52.0)
FORMAL_FR: String (35.0)
NAME_CIAWF: String (45.0)
NOTE_ADM0: String (22.0)
NOTE_BRK: String (63.0)
NAME_SORT: String (36.0)
NAME_ALT: String (19.0)
MAPCOLOR7: Integer (1.0)
MAPCOLOR8: Integer (1.0)
MAPCOLOR9: Integer (1.0)
MAPCOLOR13: Integer (3.0)
POP_EST: Integer64 (10.0)
POP_RANK: Integer (2.0)
GDP_MD_EST: Real (11.2)
POP_YEAR: Integer (4.0)
LASTCENSUS: Integer (4.0)
GDP_YEAR: Integer (4.0)
ECONOMY: String (26.0)
INCOME_GRP: String (23.0)
WIKIPEDIA: Integer (3.0)
FIPS_10_: String (3.0)
ISO_A2: String (3.0)
ISO_A3: String (3.0)
ISO_A3_EH: String (3.0)
ISO_N3: String (3.0)
UN_A3: String (4.0)
WB_A2: String (3.0)
WB_A3: String (3.0)
WOE_ID: Integer (8.0)
WOE_ID_EH: Integer (8.0)
WOE_NOTE: String (167.0)
ADM0_A3_IS: String (3.0)
ADM0_A3_US: String (3.0)
ADM0_A3_UN: Integer (3.0)
ADM0_A3_WB: Integer (3.0)
CONTINENT: String (23.0)
REGION_UN: String (23.0)
SUBREGION: String (25.0)
REGION_WB: String (26.0)
NAME_LEN: Integer (2.0)
LONG_LEN: Integer (2.0)
ABBREV_LEN: Integer (2.0)
TINY: Integer (3.0)
HOMEPART: Integer (3.0)
MIN_ZOOM: Real (3.1)
MIN_LABEL: Real (3.1)
MAX_LABEL: Real (4.1)
NE_ID: Integer64 (10.0)
WIKIDATAID: String (8.0)
NAME_AR: String (72.0)
NAME_BN: String (148.0)
NAME_DE: String (46.0)
NAME_EN: String (44.0)
NAME_ES: String (44.0)
NAME_FR: String (54.0)
NAME_EL: String (88.0)
NAME_HI: String (126.0)
NAME_HU: String (52.0)
NAME_ID: String (46.0)
NAME_IT: String (48.0)
NAME_JA: String (63.0)
NAME_KO: String (47.0)
NAME_NL: String (49.0)
NAME_PL: String (47.0)
NAME_PT: String (43.0)
NAME_RU: String (86.0)
NAME_SV: String (57.0)
NAME_TR: String (42.0)
NAME_VI: String (56.0)
NAME_ZH: String (36.0)

ogrinfo -so命令列出了shape file的总结信息，可以看到数据的主要信息（类似于GeoPandas的DataFrame）。类似地，也可以用ogr将shapefile转换为GeoJSON文件，可以用如下命令：

1
!ogr2ogr -f "GeoJSON" "data/ne_10m_admin_0_countries/ogr_data.json" "data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp"

此外，OGR也可以用来读写KML文件。

Pythonic做法

前面实际上使用的是OGR的命令行，也可以直接用Python的代码来实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
from osgeo import ogr

source = ogr.Open('data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp')
layer = source.GetLayer()
schema = []
ldefn = layer.GetLayerDefn()
for i in range(ldefn.GetFieldCount()):
    fdefn = ldefn.GetFieldDefn(i)
    schema.append(fdefn.name)


print(schema)
print(len(schema))
['featurecla', 'scalerank', 'LABELRANK', 'SOVEREIGNT', 'SOV_A3', 'ADM0_DIF', 'LEVEL', 'TYPE', 'ADMIN', 'ADM0_A3', 'GEOU_DIF', 'GEOUNIT', 'GU_A3', 'SU_DIF', 'SUBUNIT', 'SU_A3', 'BRK_DIFF', 'NAME', 'NAME_LONG', 'BRK_A3', 'BRK_NAME', 'BRK_GROUP', 'ABBREV', 'POSTAL', 'FORMAL_EN', 'FORMAL_FR', 'NAME_CIAWF', 'NOTE_ADM0', 'NOTE_BRK', 'NAME_SORT', 'NAME_ALT', 'MAPCOLOR7', 'MAPCOLOR8', 'MAPCOLOR9', 'MAPCOLOR13', 'POP_EST', 'POP_RANK', 'GDP_MD_EST', 'POP_YEAR', 'LASTCENSUS', 'GDP_YEAR', 'ECONOMY', 'INCOME_GRP', 'WIKIPEDIA', 'FIPS_10_', 'ISO_A2', 'ISO_A3', 'ISO_A3_EH', 'ISO_N3', 'UN_A3', 'WB_A2', 'WB_A3', 'WOE_ID', 'WOE_ID_EH', 'WOE_NOTE', 'ADM0_A3_IS', 'ADM0_A3_US', 'ADM0_A3_UN', 'ADM0_A3_WB', 'CONTINENT', 'REGION_UN', 'SUBREGION', 'REGION_WB', 'NAME_LEN', 'LONG_LEN', 'ABBREV_LEN', 'TINY', 'HOMEPART', 'MIN_ZOOM', 'MIN_LABEL', 'MAX_LABEL', 'NE_ID', 'WIKIDATAID', 'NAME_AR', 'NAME_BN', 'NAME_DE', 'NAME_EN', 'NAME_ES', 'NAME_FR', 'NAME_EL', 'NAME_HI', 'NAME_HU', 'NAME_ID', 'NAME_IT', 'NAME_JA', 'NAME_KO', 'NAME_NL', 'NAME_PL', 'NAME_PT', 'NAME_RU', 'NAME_SV', 'NAME_TR', 'NAME_VI', 'NAME_ZH']
94

与GeoPandas相比，OGR的操作更加繁琐一些，需要首先获得图层，然后迭代得到每个图层的属性名。

对于每一个对象，我们可以通过属性名，查看该对象的相应属性，例如我们以ADMIN为例，可以查看地图上的国家/地区名。此外，GetFeatureCount()则统计了shp文件中对象的个数。

1
2
3
4
for i in range(5):
    feature = layer[i]
    print(feature.GetField("ADMIN"))
print(layer.GetFeatureCount())
Indonesia
Malaysia
Chile
Bolivia
Peru
255

如前所述，CRS是空间数据中非常重要的信息，我们可以通过两种方式获得这一信息——直接读取图层和从图层的几何图形中，具体代码如下：

1
2
3
4
5
6
7
8
9
# Option 1: from Layer
spatialRef = layer.GetSpatialRef()
print(spatialRef)
# Option 2: from Geometry
layer = source.GetLayer()
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef2 = geom.GetSpatialReference()
print(spatialRef2)
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

如上所示，GetNextFeature()获得图形对象，GetGeometryRef()获得“句柄”，可以进一步对每个图形进行操作和查看，例如查看对象的中心点：

1
print(geom.Centroid())
POINT (117.270433339165 -2.22296100251739)

读写栅格数据

Rasterio库方法

导入Rasterio库后，打开一个GeoTIFF文件，然后可以用一些简单的命令查看数据。

1
import rasterio
1
dataset = rasterio.open('data/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif')
1
2
3
4
5
6
print(f'The number of image bans is {dataset.count}.')
print(f'The number of columns is {dataset.width}.')
print(f'The number of rows is {dataset.height}.')
print(f'The spatial bounding box is {dataset.bounds}.')
print(f'The CRS of the dataset is {dataset.crs}.')

The number of image bans is 3.
The number of columns is 21600.
The number of rows is 10800.
The spatial bounding box is BoundingBox(left=-180.0, bottom=-90.000000000036, right=180.00000000007202, top=90.00000000000001).
The CRS of the dataset is EPSG:4326.

对于栅格数据的不同光谱，可以用read()返回ndarray类型的数据：

1
2
band1 = dataset.read(1)
band1
array([[110, 110, 108, ..., 107, 106, 106],
       [110, 110, 108, ..., 107, 106, 106],
       [107, 107, 105, ..., 106, 107, 107],
       ...,
       [230, 230, 231, ..., 231, 230, 230],
       [231, 231, 231, ..., 231, 230, 231],
       [230, 230, 231, ..., 229, 230, 230]], dtype=uint8)

使用matplotlib将栅格数据可视化：

1
2
import matplotlib.pyplot as plt
%matplotlib inline
1
2
3
plt.figure(figsize = (15, 12))
plt.imshow(dataset.read(1))
plt.show()

output_61_0

GDAL库方法

类似矢量数据的读取方法，用gdalinfo读写栅格数据：

1
!gdalinfo "data/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif"
Driver: GTiff/GeoTIFF
Files: data/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif
Size is 21600, 10800
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000014)
Pixel Size = (0.016666666666670,-0.016666666666670)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2014:10:18 12:06:24
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CC 2014 (Macintosh)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (   0.0000000,  -0.0000000) (  0d 0' 0.00"E,  0d 0' 0.00"S)
Band 1 Block=21600x1 Type=Byte, ColorInterp=Red
Band 2 Block=21600x1 Type=Byte, ColorInterp=Green
Band 3 Block=21600x1 Type=Byte, ColorInterp=Blue

用gdal_translate命令将GeoTFIFF文件转为JPEG格式：

1
!gdal_translate -of JPEG "data/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif" output.jpg
Input file size is 21600, 10800
0...10...20...30...40...50...60...70...80...90...100 - done.

1
2
3
4
5
6
from PIL import Image
Image.MAX_IMAGE_PIXELS = 1000000000
img = Image.open('output.jpg')
plt.figure(figsize = (15, 12))
plt.imshow(img)
plt.show()

output_67_0

GDAL打开GeoPackage

GeoPackage既可以是矢量数据，也可以是栅格数据。首先用gdalinfo查看GeoPackage数据的信息：

1
!gdalinfo "data/gdal_sample_v1.2_no_extensions.gpkg"
Driver: GPKG/GeoPackage
Files: data/gdal_sample_v1.2_no_extensions.gpkg
Size is 512, 512
Coordinate System is
Subdatasets:
  SUBDATASET_1_NAME=GPKG:data/gdal_sample_v1.2_no_extensions.gpkg:byte_png
  SUBDATASET_1_DESC=byte_png - byte_png
  SUBDATASET_2_NAME=GPKG:data/gdal_sample_v1.2_no_extensions.gpkg:byte_jpeg
  SUBDATASET_2_DESC=byte_jpeg - byte_jpeg
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

GDAL下载数据

GDAL Web地图服务（WMS）驱动程序允许与在线Web地图服务进行交互。 您可以使用它直接从命令提示符（或在本例中为Jupyter Notebook）下载各种地理空间数据集，子集或有关可用数据集的信息，而无需使用浏览器导航到网站并手动下载数据。 有许多不同的选项，因此请参阅联机文档以获取更多信息。 以下示例需要GDAL 2.0或更高版本。 以下命令使用ArcGIS MapServer的REpresentational State Transfer（REST）定义的URL，并返回有关所请求的图像服务的信息，例如波段数，波段名称，CRS，角坐标等：

1
!gdalinfo "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json&pretty=true"
Driver: WMS/OGC Web Map Service
Files: none associated
Size is 1073741824, 1073741824
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+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 +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-20037508.342787001281977,20037508.342787001281977)
Pixel Size = (0.037322767717361,-0.037322767717361)
Metadata:
  CACHE_PATH=./gdalwmscache\1e1be7907a3b947a2d7ff2cd40868eb5
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-20037508.343,20037508.343) (180d 0' 0.00"W, 85d 3' 4.06"N)
Lower Left  (-20037508.343,-20037508.343) (180d 0' 0.00"W, 85d 3' 4.06"S)
Upper Right (20037508.343,20037508.343) (180d 0' 0.00"E, 85d 3' 4.06"N)
Lower Right (20037508.343,-20037508.343) (180d 0' 0.00"E, 85d 3' 4.06"S)
Center      (  -0.0000030,   0.0000030) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 536870912x536870912, 268435456x268435456, 134217728x134217728, 67108864x67108864, 33554432x33554432, 16777216x16777216, 8388608x8388608, 4194304x4194304, 2097152x2097152, 1048576x1048576, 524288x524288, 262144x262144, 131072x131072, 65536x65536, 32768x32768, 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512, 256x256
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 536870912x536870912, 268435456x268435456, 134217728x134217728, 67108864x67108864, 33554432x33554432, 16777216x16777216, 8388608x8388608, 4194304x4194304, 2097152x2097152, 1048576x1048576, 524288x524288, 262144x262144, 131072x131072, 65536x65536, 32768x32768, 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512, 256x256
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 536870912x536870912, 268435456x268435456, 134217728x134217728, 67108864x67108864, 33554432x33554432, 16777216x16777216, 8388608x8388608, 4194304x4194304, 2097152x2097152, 1048576x1048576, 524288x524288, 262144x262144, 131072x131072, 65536x65536, 32768x32768, 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512, 256x256

小结

本文简要介绍了GIS中主要的数据类型，解释了矢量和栅格数据的差异，然后用GeoPandas、OGR、GDAL和Rasterio等库进行了地理空间数据的读写示例。接下来，我们将利用GGR、Shapely和GeoPandas等库进行地理空间分析和处理。