ほぼコピペですが、頑張ります。
https://www.kaggle.com/code/alexisbcook/coordinate-reference-systems
導入
このコースで作成する地図は、地球の表面を二次元で表現しています。しかし、実際の世界は三次元の地球です。そのため、地球を平面として表示するためにマップ投影と呼ばれる方法を使用する必要があります。
マップ投影は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)
座標参照系は、European Petroleum Survey Group(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の両方を使用する必要がありました:
1. 緯度と経度の座標を含む列を持つDataFrameを作成します。
2. それをGeoDataFrameに変換するために、`gpd.GeoDataFrame()` を使用します。
3. `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)
`to_crs()` メソッドは「geometry」列のみを変更し、他のすべての列はそのままにします。
# The "Latitude" and "Longitude" columns are unchanged
facilities.to_crs(epsg=32630).head()
GeoPandasでEPSGコードが利用できない場合、CRSを「proj4文字列」として知られるもので変更することができます。たとえば、緯度/経度座標に変換するためのproj4文字列は以下の通りです:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
In [6]:
# Change the CRS to EPSG 4326
regions.to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs").head()
幾何学的オブジェクトの属性
最初のチュートリアルで学んだように、任意のGeoDataFrameにおいて、「geometry」列のタイプは表示しようとしているものに依存します。たとえば、次のように使用できます:
- 地震の震源地のためにPoint
- 道路のためにLineString
- 国境を示すためにPolygon
これらの3つの種類の幾何学的オブジェクトには、データセットを迅速に分析するために使用できる組み込み属性があります。たとえば、Pointのx属性とy属性を使用して、Pointのx座標とy座標を取得できます。
# Get the x-coordinate of each point
facilities.geometry.head().x
Out[7]:
0 -1.96317
1 -1.58592
2 -1.34982
3 -1.61098
4 -1.61098
dtype: float64
また、LineStringの長さは`length`属性から取得できます。
また、Polygonの面積は`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()
Area of Ghana: 239584.5760055668 square kilometers
CRS: PROJCS["WGS_1984_UTM_Zone_30N",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]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
上記のコードセルでは、地域のGeoDataFrameのCRSがEPSG 32630(「メルカトル」投影)に設定されているため、面積の計算は「アフリカアルバーズ等積円錐」のような等積投影を使用した場合よりもわずかに正確さに欠けるかもしれません。
しかし、これにより、ガーナの面積は約239,585平方キロメートルとなり、正確な答えからそれほど大きくはずれていません。