基于Python的地理空间分析(三):GIS数据类型介绍

前面已经介绍了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
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.png

查看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.png

如果想将读取的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
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.png

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.png

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等库进行地理空间分析和处理。

张da统帅 wechat
扫码订阅我的微信公众号