いきなり長文です。
A coordinate reference systems (CRS) is a crucial piece of metadata for any geospatial data set. Without a CRS, the geometries would simply be a collection of coordinates in an arbitrary space. Only the CRS allows GIS software, including the Python packages we use in this course, to relate these coordinates to a place on Earth (or other approximately spherical objects or planets).
ざっくりというと、座標というデータを特定の場所に関連付けるためには、CRSが必要です。
Often conflated with coordinate reference systems, and definitely closely related, are map projections. Map projections, also called projected coordinate systems, are mathematical models that allow us to transfer coordinates on the surface of our three-dimensional Earth into coordinates in a planar surface, such as a flat, two-dimensional map. In contrast to projected coordinate systems, geographic coordinate systems simply directly use latitude and longitude, i.e. the degrees along the horizontal and vertical great circles of a sphere approximating the Earth, as the x and y coordinates in a planar map. Finally, there are both projected and geographic coordinate systems that make use of more complex ellipsoids than a simple sphere to better approximate the ‘potato-shaped’ reality of our planet. The full CRS information needed to accurately relate geospatial information to a place on Earth includes both (projected/geographic) coordinate system and ellipsoid.
CRSとは別物ですが、関連するものとして、map projections、投影座標系があります。
三次元である地球の表面の座標を**平面（二次元）**の座標に転送できる数学的モデルです。
しばしば座標参照システムと混同され、間違いなく密接に関連している地図投影です。地図の投影は、投影座標系とも呼ばれ、3次元の地球の表面の座標を平面の平面の座標に転送できる数学的モデルです。
それは別に、地理座標系というものがあります。latitudeとlongitude、いわるゆ緯度・経度です。
投影座標系・地理座標系いずれも、地球そのものに近づけるため、単純な球体ではなく、より複雑な楕円体を用います。
The CRS in different spatial datasets differ fairly often, as different coordinate systems are optimised for certain regions and purposes. No coordinate system can be perfectly accurate around the globe, and the transformation from three- to two-dimensional coordinates can not be accurate in angles, distances, and areas simultaneously.
座標系は特定の領域と目的に最適化されていいます。世界中どこでも正確な座標系は無く、また、三次元座標から二次元座標への変換において、角度・距離・面積のすべてを正確に満たす座標系もありません。
Consequently, it is a common GIS task to transform (or reproject) a data set from one references system into another, for instance, to make two layers interoperatable. Comparing two data sets that have different CRS would inevitably produce wrong results; for example, finding points contained within a polygon cannot work, if the the points have geographic coordinates (in degrees), and the polygon is in the national Finnish reference system (in meters).
なので、データセットを組み合わせるときは、その座標系のつじつま合わせが必要になります。
Choosing an appropriate projection for your map is not always straightforward. It depends on what you actually want to represent in your map, and what your data’s spatial scale, resolution and extent are. In fact, there is not a single ‘perfect projection’; each has strengths and weaknesses, and you should choose a projection that fits best for each map. In fact, the projection you choose might even tell something about you:
ただ、投影座標系を選ぶことは、それほど簡単ではありません。それぞれの投影座標系は長所・短所があるからでです。
Handling coordinate reference systems in Geopandas
Once you have figured out which map projection to use, handling coordinate reference systems, fortunately, is fairly easy in Geopandas. The library pyproj provides additional information about a CRS and can assist with more tricky tasks, such as guessing the unknown CRS of a data set.
In this section, we will learn how to retrieve the coordinate reference system information of a data set, and how to re-project the data into another CRS.
GeopandasでのCRSの操作は簡単です。pyproj ライブラリを使用します。
では、データセットのCRSの取得と、別のCRSへの変換（再投影）をやってみましょう。
Displaying the CRS of a data set
We will start by loading a data set of EU countries that has been downloaded from the Geographic Information System of the Commission (GISCO), a unit within Eurostat that manages geospatial data for the European Commission.
EU各国のデータセットしようします。ファイルはここにあります。
https://github.com/Automating-GIS-processes/site/tree/main/docs/lessons/lesson-2/data/eu_countries
Let’s check the data set’s CRS:
では、CRSを見てみましょう。
import pathlib
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"
import geopandas
eu_countries = geopandas.read_file(
DATA_DIRECTORY / "eu_countries" / "eu_countries_2022.gpkg"
)
eu_countries.crs
What we see here is, in fact, a pyproj.CRS object.
The EPSG code (European Petroleum Survey Group) is global standard for identifying coordinate reference systems. The number refers to an entry in the EPSG Geodetic Parameter Dataset, a collection of coordinate reference systems coordinate transformations ranging from global to national, regional, and local scope.
Our GeoDataFrame’s EPSG code is 4326. This is number to remember, as you will come across it a lot in the geospatial world: It refers to a geographic coordinate system using the WGS-84 reference ellipsoid. This is the most commonly used coordinate reference system in the world. It’s what we refer to when we speak of longitude and latitude.
You can find information about reference systems and lists of commonly known CRS from many online resources, for example: - spatialreference.org - proj4.org - mapref.org
参照しているのは、pyproj.CRSのオブジェクトです。EPSG（European Petroleum Survey Group）コードは世界基準のCRSです。EPSG:4326は頻出で、WGS-84参照楕円体による地理的座標系です。
Indeed, the coordinate values of our data set seem to be in an appropriate value range.
EPSG:4326では、緯度は-180度から180度、緯度は-90度から90度であり、実際のデータもそうなっていることが確認できます。
Reprojecting a GeoDataFrame
A geographic coordinate system, EPSG:4326, is not particularly well-suited for showing the countries of the European Union. Distortion is high. Rather, we could use a Lambert Azimuthal Equal-Area projection, such as
EPSG:3035 <https://spatialreference.org/ref/epsg/etrs89-etrs-laea/>__, the map projection officially recommended by the European Commission.
Transforming data from one reference system to another is a very simple task in geopandas. In fact, all you have to to is use the to_crs() method of a GeoDataFrame, supplying a new CRS in a wide range of possible formats. The easiest is to use an EPSG code:
EPSG:4326でEUの各国を表示すると歪みが大きい。欧州委員会が公式に推奨しているEPSG:3035に変換します。
GeoDataFrameのto_crs()メソッドの引数にEPSGコードを設定します。
eu_countries_EPSG3035 = eu_countries.to_crs("EPSG:3035")
eu_countries_EPSG3035.geometry.head()
Let’s check how the coordinate values have changed:
And here we go, the coordinate values in the geometries have changed! Congratulations on carrying out your very first geopandas coordinate transformation!
To better grasp what exactly we have just done, it is a good idea to explore our data visually. Let’s plot our data set both before and after the coordinate transformation. We will use matplotlib’s subplots feature that we got to know in Geo-Python lesson 7.
座標変換の前後のデータセットをプロットしましょう。それがわかりやすいです。
import matplotlib.pyplot
# Prepare sub plots that are next to each other
figure, (axis1, axis2) = matplotlib.pyplot.subplots(nrows=1, ncols=2)
# Plot the original (WGS84, EPSG:4326) data set
eu_countries.plot(ax=axis1)
axis1.set_title("WGS84")
axis1.set_aspect(1)
# Plot the reprojected (EPSG:3035) data set
eu_countries_EPSG3035.plot(ax=axis2)
axis2.set_title("ETRS-LAEA")
axis2.set_aspect(1)
matplotlib.pyplot.tight_layout()
Indeed, the maps look quite different, and the re-projected data set distorts the European countries less, especially in the Northern part of the continent.
明らかに違いあります。EPSG:3035（右側）はヨーロッパ大陸、特に北部の歪みが少ないです。
Let’s still save the reprojected data set in a file so we can use it later. Note that, even though modern file formats save the CRS reliably, it is a good idea to use a descriptive file name that includes the reference system information.
EPSG:3035のデータをファイルに保存しましょう。あとで使うこともできますので。ファイル名でCRSがわかるようにしておくと良いでしょう。
eu_countries_EPSG3035.to_file(
DATA_DIRECTORY / "eu_countries" / "eu_countries_EPSG3035.gpkg"
)
Handling different CRS formats
There are different ways to store and represent CRS information. The more commonly used formats include PROJ strings, EPSG codes, Well-Known-Text (WKT) and JSON. You will likely encounter some or all of these when working with spatial data obtained from different sources. Being able to convert CRS information from one format to another is needed every now and then, hence, it is useful to know a few tricks how to do this.
We’ve already briefly mentioned that geopandas uses the pyproj library to handle reference systems. We can use the same module to parse and convert CRS information in different formats.
CRSについてのtipsです。
Overview
Below we print different representations of the CRS of the data set of EU countries we used before:
いろいろな形式でCRSを表示してみましょう。
import pyproj
crs = pyproj.CRS(eu_countries.crs)
print(f"CRS as a proj4 string: {crs.to_proj4()}\n")
print(f"CRS in WKT format: {crs.to_wkt()}\n")
print(f"EPSG code of the CRS: {crs.to_epsg()}\n")
Use pyproj to find detailed information about a CRS
A pyproj.CRS object can also be initialised manually, for instance, using an EPSG code or a Proj4-string. It can then provide detailed information on the parameters of the reference system, as well as suggested areas of use. We can, for example, create a CRS object for the EPSG:3035 map projection we used above:
pyprojの機能で、CRSの詳細を見てみましょう。
crs = pyproj.CRS("EPSG:4326")
crs
Immediately, we see that a pyproj.CRS object contains rich information about the reference system: a name, an area of use (including bounds), a datum and an ellipsoid.
pyproj.CRSオブジェクトには、名前、使用領域、データム、楕円体など、様々な情報が含まれています。
この情報は、個別に抽出することもできます。
This information can also be extraced individually:
Global map projections
Finally, it’s time to play around with some map projections. For this, you will find a global data set of country polygons in the data directory. It was downloaded from naturalearthdata.com, a fantastic resource for cartographer-grade geodata.
最後に、CRSの変換（再投影）をいくつかやってみましょう。地球全体のデータセットが必要です。ファイルは以下にあります。
データをそのまま表示しました。CRSはWGS-84でした。
world_countries = geopandas.read_file(
DATA_DIRECTORY / "world_countries" / "ne_110m_admin_0_countries.zip"
)
world_countries.plot()
matplotlib.pyplot.title(world_countries.crs.name)
CRSをEPSG:3857、Webメルカトル（Pseudo-Mercator）に変換して表示します。
Web メルカトル（EPSG:3857）について整理する
南極が愉快なことになっていますね。
# web mercator
world_countries_EPSG3857 = world_countries.to_crs("EPSG:3857")
world_countries_EPSG3857.plot()
matplotlib.pyplot.title(world_countries_EPSG3857.crs.name)
# remove axis decorations
matplotlib.pyplot.axis("off")
CRSをエケルト第4図法（Eckert IV）、EPSG:53012に変換して表示します。to_crsメソッドには、そのCRSの構成要素の文字列を指定することもできます。
エケルト図法 (第 4 図法)
外側の子午線は半円のため、図法は丸みを帯びた形状となり、外側の子午線と極線が交わる角は滑らかになります。 この図法は、正確な面積を必要とする主題世界地図やその他の世界地図によく使用されます。
# Eckert-IV (https://spatialreference.org/ref/esri/54012/)
ECKERT_IV = "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
world_countries_eckert_iv = world_countries.to_crs(ECKERT_IV)
world_countries_eckert_iv.plot()
matplotlib.pyplot.title("Eckert Ⅳ")
matplotlib.pyplot.axis("off")
最後に、独自のCRSを指定して表示します。
http://www.statsmapsnpix.com/2019/09/globe-projections-and-insets-in-qgis.html の内容をもとに、「フィンランドを中心になるように地球を真上から見たような」投影をしているようです。
# An orthographic projection, centered in Finland!
# (http://www.statsmapsnpix.com/2019/09/globe-projections-and-insets-in-qgis.html)
world_countries_ortho = world_countries.to_crs(
"+proj=ortho +lat_0=60.00 +lon_0=23.0000 +x_0=0 +y_0=0 "
"+a=6370997 +b=6370997 +units=m +no_defs"
)
world_countries_ortho.plot()
matplotlib.pyplot.axis("off")
