kaggle : Geospatial Analysis ②

본 포스팅은 kaggle : Geospatial Analysis을 수료하고 정리한 글입니다.

---

2. Coordinate Reference Systems

---

2.1 Introduction

---

지구의 표면을 2차원으로 묘사해서 지도에 나타내지만, 실제 지구는 3차원 구체이다. 그래서 map projection 방법을 사용하여 평면 표면으로 렌더링 해야한다. 한 마디로 3차원 지구타원체를 2차원 평면의 네모난 지도에 펼쳐야 한다!

렌더링이란? 3차원 공간에 객체(Object)를 2차원 화면의 하나의 장면(scene)에 바꾸어 표현하는 것을 의미한다.

각각의 Map projection 방법은 지구 표면을 알 수 없는(?) 방식으로 왜곡하지만, 동시에 유용한 속성을 유지한다. 예를 들면

• 면적 보존 투영(the equal-area projection) : 면적을 보존한다. 국가나 도시의 면적을 계산하려는 경우에 효과적이다.
• 등거리 투영(the equidistant projection) : 거리를 보존하다. 비행 거리를 계산하는데 쓸모가 있다.

좌표 참조 시스템(CRS)을 사용하여 투영된 점들이 지구상의 실제 위치와 어떻게 부합하는지 볼 것이다. 그리고 GeoPandas에서 좌표 참조 시스템을 사용하는 방법에 대해서 알아보자.

2.2 Setting the CRS

---

shapefile에서 GeoDataFrame을 생성하면, 좌표 참조 시스템이 import된다.

# 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)

epsg:32630

위의 epsg:32630는 각도를 보존해 해상 항법에 유용하지만 면적을 약간 왜곡시킨다는 단점이 있다.

EPSG는 European Petroleum Survey Group의 약자로, 지리 정보 시스템(GIS)에서 사용되는 좌표 참조 시스템(CRS)을 식별하기 위한 코드 체계이다. 이 체계는 전 세계적으로 사용되며, 각 CRS에는 고유한 EPSG 코드가 할당된다. 네이버지도는 'EPSG:5179' 이며, 카카오맵은 'EPSG:5181'이다. 코드별로 원점도 다르고 특성도 다르다.

만약 shp파일이 아닌 csv파일로 GeoDataFrame을 생성할 때는 CRS를 설정해야 한다.

아래 코드의 EPSG:4326은 지구의 위도와 경도를 사용한다. GPS 데이터와 같은 위치 데이터를 처리하는데 효과적이다.

# 데이터 로딩
facilities_df = pd.read_csv("../input/geospatial-learn-course-data/ghana/ghana/health_facilities.csv")

# DataFrame 에서 GeoDataFrame로 변환
facilities = gpd.GeoDataFrame(facilities_df, geometry=gpd.points_from_xy(facilities_df.Longitude, facilities_df.Latitude))

# CRS -> EPSG 4326
facilities.crs = {'init': 'epsg:4326'}

# 상위 5개 행 출력
facilities.head()

	Region	District	FacilityName	Type	Town	Ownership	Latitude	Longitude	geometry
0	Ashanti	Offinso North	A.M.E Zion Clinic	Clinic	Afrancho	CHAG	7.40801	-1.96317	POINT (-1.96317 7.40801)
1	Ashanti	Bekwai Municipal	Abenkyiman Clinic	Clinic	Anwiankwanta	Private	6.46312	-1.58592	POINT (-1.58592 6.46312)
2	Ashanti	Adansi North	Aboabo Health Centre	Health Centre	Aboabo No 2	Government	6.22393	-1.34982	POINT (-1.34982 6.22393)
3	Ashanti	Afigya-Kwabre	Aboabogya Health Centre	Health Centre	Aboabogya	Government	6.84177	-1.61098	POINT (-1.61098 6.84177)
4	Ashanti	Kwabre	Aboaso Health Centre	Health Centre	Aboaso	Government	6.84177	-1.61098	POINT (-1.61098 6.84177)

 

코드 셀을 해석해보자.

• 먼저 cvs파일을 로딩 후 DataFrame에서 GeoDataFrame으로 변환했다.
• csv파일로 GeoDataFrame을 생성했기에 CRS를 설정해주었다.
• gdp.points_from_xy()은 위도와 경도 열에서 Point를 생성한다.

2.3 Re-projecting

---

재투영(Re-projecting) 은 좌표 참조 시스템을 변경하는 과정을 말한다. GeoPandas에서는 to_crs() 메소드를 사용하여 수행한다.

두 개 이상의 GeoDataFrame을 그릴 때, 모든 GeoDataFrame이 같은 CRS를 사용하는 것이 중요하다. 다른 좌표를 가진 GeoDataFrame를 지도 위에 표시할 경우, 데이터가 지도 상에 잘못된 위치에 표시될 수도 있기 때문이다.

아래 코드는 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()

	Region	District	FacilityName	Type	Town	Ownership	Latitude	Longitude	geometry
0	Ashanti	Offinso North	A.M.E Zion Clinic	Clinic	Afrancho	CHAG	7.40801	-1.96317	POINT (614422.662 818986.851)
1	Ashanti	Bekwai Municipal	Abenkyiman Clinic	Clinic	Anwiankwanta	Private	6.46312	-1.58592	POINT (656373.863 714616.547)
2	Ashanti	Adansi North	Aboabo Health Centre	Health Centre	Aboabo No 2	Government	6.22393	-1.34982	POINT (682573.395 688243.477)
3	Ashanti	Afigya-Kwabre	Aboabogya Health Centre	Health Centre	Aboabogya	Government	6.84177	-1.61098	POINT (653484.490 756478.812)
4	Ashanti	Kwabre	Aboaso Health Centre	Health Centre	Aboaso	Government	6.84177	-1.61098	POINT (653484.490 756478.812)

EPSG 코드가 GeoPandas에 없을 경우에는, 좌표 참조 시스템의 "pro4 string"을 사용하여 CRS를 변경하면 된다. 예를 들어 epsg:4362을 변환하는 proj4 문자열은 다음과 같다.

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

	Region	geometry
0	Ashanti	POLYGON ((-1.30985 7.62302, -1.30786 7.62198, ...
1	Brong Ahafo	POLYGON ((-2.54567 8.76089, -2.54473 8.76071, ...
2	Central	POLYGON ((-2.06723 6.29473, -2.06658 6.29420, ...
3	Eastern	POLYGON ((-0.21751 7.21009, -0.21747 7.20993, ...
4	Greater Accra	POLYGON ((0.23456 6.10986, 0.23484 6.10974, 0....

2.4 Attributes of geometric objects

---

geometry 열의 유형에 따라 해당 데이터를 나타내는 방식 또한 달라지게 된다. 예를 들어

• 지진의 진앙 지점을 나타내는 경우는 Point
• 도로를 나타내는 경우는 LineString
• 국가 경계를 나타내는 경우는 Polygon

이 사용된다.

이러한 지리 객체 유형에는 내장된 속성이 있어 빠르게 데이터를 분석할 수 있다.

# Get the x-coordinate of each point
facilities.geometry.head().x

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

이렇게 x속성을 가져와 x좌표로 사용할 수 있다.

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: epsg:32630

	Region	geometry	AREA
0	Ashanti	POLYGON ((686446.075 842986.894, 686666.193 84...	24379.017777
1	Brong Ahafo	POLYGON ((549970.457 968447.094, 550073.003 96...	40098.168231
2	Central	POLYGON ((603176.584 695877.238, 603248.424 69...	9665.626760
3	Eastern	POLYGON ((807307.254 797910.553, 807311.908 79...	18987.625847
4	Greater Accra	POLYGON ((858081.638 676424.913, 858113.115 67...	3706.511145

 

등면적 투영법이 아니라 면적 계산에 약간의 오차는 있지만 실제 면적과 크게 차이나지 않아 충분히 분석에 활용가능하다.