时间戳

首页 > 技术文章 > kaggle地理空间分析-coordinate reference system(坐标参考系统)<ii>

kaggle地理空间分析-coordinate reference system(坐标参考系统)<ii>

hy1231 2020-08-08 18:56 原文

介绍

您在本课程中创建的地图以二维方式描绘了地球表面。但是，正如您所知，世界实际上是三维地球。因此，我们必须使用一种称为地图投影的方法将其渲染为平坦表面。

地图投影不能100％准确。每个投影都会以某种方式扭曲地球表面，同时保留一些有用的属性。例如，

  • 等面积投影（例如“兰伯特圆柱等分面积”或“非洲阿尔伯斯等分圆锥”）可保留面积。例如，如果您想计算一个国家或城市的面积，这是一个不错的选择。

  • 等距投影（例如“方位角等距投影”）保留距离。这将是计算飞行距离的好选择。

我们使用坐标参考系统（CRS）来显示投影点如何对应于地球上的真实位置。在本教程中，您将了解有关坐标参考系统的更多信息，以及如何在GeoPandas中使用它们。

import geopandas as gpd
import pandas as pd

设置CRS

当我们从shapefile创建GeoDataFrame时，已经为我们导入了CRS。

# Load a GeoDataFrame containing regions in Ghana
regions = gpd.read_file("../input/geospatial-learn-course-data/ghana/ghana/Regions/Map_of_Regions_in_Ghana.shp")
print(regions.crs)

{'init': 'epsg:32630'}

您如何解释？

坐标参考系统由欧洲石油测量集团（EPSG）代码参考。

此GeoDataFrame使用EPSG 32630，通常更称为“墨卡托”投影。 该投影保留了角度（使其对于海上航行很有用）并且使区域略微变形。

但是，从CSV文件创建GeoDataFrame时，必须设置CRS。 EPSG 4326对应于纬度和经度的坐标。

# Create a DataFrame with health facilities in Ghana
facilities_df = pd.read_csv("../input/geospatial-learn-course-data/ghana/ghana/health_facilities.csv")

# Convert the DataFrame to a GeoDataFrame
facilities = gpd.GeoDataFrame(facilities_df, geometry=gpd.points_from_xy(facilities_df.Longitude, facilities_df.Latitude))

# Set the coordinate reference system (CRS) to EPSG 4326
facilities.crs = {'init': 'epsg:4326'}

# View the first five rows of the GeoDataFrame
facilities.head()

在上面的代码单元中，要从CSV文件创建GeoDataFrame，我们需要同时使用Pandas和GeoPandas：

  • 我们首先创建一个DataFrame，其中包含具有经度和纬度坐标的列。
  • 要将其转换为GeoDataFrame，我们使用gpd.GeoDataFrame（）。
  • gpd.points_from_xy（）函数从纬度和经度列创建Point对象。

重新投影

重新投影是指更改CRS的过程。 这是在GeoPandas中使用to_crs（）方法完成的。

绘制多个GeoDataFrame时，重要的是它们都使用相同的CRS。 在下面的代码单元中，我们在绘制之前更改设施GeoDataFrame的CRS以匹配区域的CRS。

# Create a map
ax = regions.plot(figsize=(8,8), color='whitesmoke', linestyle=':', edgecolor='black')
facilities.to_crs(epsg=32630).plot(markersize=1, ax=ax)

<matplotlib.axes._subplots.AxesSubplot at 0x7faa424e6748>


to_crs（）方法仅修改“ geometry”列：所有其他列均保持不变。

# The "Latitude" and "Longitude" columns are unchanged
facilities.to_crs(epsg=32630).head()

如果在GeoPandas中没有EPSG代码，我们可以使用CRS的“ proj4字符串”来更改CRS。 例如，要转换为纬度/经度坐标的proj4字符串如下：

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

# Change the CRS to EPSG 4326
regions.to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs").head()

几何对象的属性

正如您在第一个教程中所了解的那样，对于任意的GeoDataFrame，“ geometry”列中的类型取决于我们要显示的内容：例如，我们可以使用：

  • 震中点
  • 街道的LineString，或
  • 显示国家边界的多边形。
    三种类型的几何对象均具有内置属性，可用于快速分析数据集。 例如，您可以分别从x和y属性获得Point的x和y坐标。
# Get the x-coordinate of each point
facilities.geometry.x.head()

out:
0   -1.96317
1   -1.58592
2   -1.34982
3   -1.61098
4   -1.61098
dtype: float64

并且，您可以从length属性中获取LineString的长度。

或者，您可以从area属性获得多边形的面积。

# Calculate the area (in square meters) of each polygon in the GeoDataFrame 
regions.loc[:, "AREA"] = regions.geometry.area / 10**6

print("Area of Ghana: {} square kilometers".format(regions.AREA.sum()))
print("CRS:", regions.crs)
regions.head()

out:
Area of Ghana: 239584.5760055668 square kilometers
CRS: {'init': 'epsg:32630'}

在上面的代码单元中，由于区域GeoDataFrame的CRS设置为EPSG 32630（“墨卡托”投影），因此与使用等面积投影（如“非洲阿尔伯斯等距圆锥形”）相比，面积计算的准确性稍差 ”。

但是，这加纳的面积约为239585平方公里，与正确答案相差不远。

推荐阅读