python으로 공간 데이터 시각화를 하는 방법은 두 가지가 있습니다.

1. matplotlib
2. plotly

과학자들은 matplotlib으로 그래프를 그리는 것이 익숙합니다.
공간 데이터 시각화에도 matplotlib 기반으로 많이 접근합니다.
plotly는 지도를 그려놓고 그 위에서 공간 데이터를 시각화합니다.
오늘 할 이야기는 matplotlib을 기반으로 공간 데이터를 시각화하는 얘기, 그 중에서도 Basemap을 사용하는 방법에 대해 알아보고자 합니다.

왜 Basemap이냐?

matplotlib기반 공간 정보 시각화 방법에는 cartopy와 basemap이 있습니다.
이때까지는 cartopy를 많이 썼습니다.
제가 matplotlib 기반으로 그래프는 많이 그려봤지만 cartopy는 그 명성은 익히 들었지만 거의 사용해 본 적이 없었습니다.
일단 무겁고, 설치 환경에서 자꾸 환경 설정에 실패했기 때문입니다.
그 대안으로 어떤 라이브러리를 사용할 까 고민하던 와중에 발견한 것이 Basemap입니다.

cartopy와 basemap의 차이점은 다음과 같습니다.
둘 다 matplotlib 기반이지만, basemap은 matplotlib에서 만든 라이브러리고 cartopy는 그 기반이 matplotlib이지만 자체 생태계를 구축했습니다.
저도 basemap이 matplotlib에서 만든 줄은 이번에 글을 쓰기 위해 자료를 찾아보던 중 알게 되었습니다.
basemap 공식 문서

어떻게 사용하는가?

GEOS Projection?

이번 포스팅에서는 GEOS Projection을 따르는 위성 영상을 시각화 하기 위해서 basemap을 사용합니다.
GEOS Projection이 뭔지에 대해서는 아래 링크를 참조하세요.
GEOS

쉽게 얘기하자면 위성 영상이 전 지구(Full Disk)를 찍은 모습을 나타내는 projection입니다.
아래 그림처럼요.

+proj=geos +h=35785831.0 +lon_0=128.2 +sweep=y

lon_0은 세로로 지구본을 얼마나 돌릴 건지 결정합니다.

설치

• 설치 링크 : https://matplotlib.org/basemap/
• pip install : https://pypi.org/project/basemap/
  • basemap-data-hires : 라이브러리 사용 시 high-resolution 자료들(해안선 등)을 지도에 쓰려면 반드시 설치해야 합니다.
  • 기본 해상도 데이터는 basemap-data로 basemap 설치할 때 같이 깔립니다.

    python -m pip install basemap
    python -m pip install basemap-data-hires

사용 예시

기본 GEOS 표출

from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt

plt.figure(figsize=(10,10))
m2 = Basemap(rsphere=(6378137.00,6356752.3142), projection='geos',lon_0=128.2,resolution='h')
m2.drawcoastlines()

원하는 영역만 표출하기

lat1, lat2, long1, long2 = (46.830415916691464,
 -16.043669898267705,
 90.29636166074948,
 172.08684795392023)

plt.figure(figsize=(10,10))
'''
llcrnr : left lower corner
urcrnr : upper right corner
resolution : basemap에서 제공하는 해안선 등의 데이터를 어떤 해상도로 표출할 지 결정. h면 high resolution을 의미
lon_0 : longitude 방향으로 중앙
lat_0 : latitude 방향으로 중앙
'''
m= Basemap(llcrnrlon=long1,llcrnrlat=lat2,urcrnrlon=long2,urcrnrlat=lat1, resolution='h',epsg=4326,lon_0=(long2-long1) / 2,lat_0=(lat1-lat2) / 2)
m.drawcoastlines()

표출 예시

NOAA에서 제공하는 GOES위성의 어떤 밴드 영상을 단순 시각화한 예시

결론

Cartopy보다 덜 무거우면서도 친숙한 인터페이스 내에 Projected map을 사용하고 싶다면 basemap 사용을 고려할 것.

그리고 기초편이 있다는 것은 심화 편도 있다는 얘기.
심화 편에서는 그렇다면 cartopy에서 그리던 것들을 어떻게 basemap 코드로 옮길 수 있는데?에 대해 쓸 예정입니다.
일하면서 틈틈히 메뉴얼을 작성하고 그 메뉴얼을 더 다듬어서 틈나는 대로 블로그 포스팅을 하다 보니 시간이 좀 걸립니다.

Reference

• geostationary projection에 대한 자세한 정보 : https://proj.org/en/9.3/operations/projections/geos.html
• basemap의 메뉴얼 중 geostationary projection에 대한 샘플 코드와 설명 : https://matplotlib.org/basemap/users/geos.html
• basemap의 여러 기본 사용법을 나열한 곳. Python Data Science Handbook 내 실려있는 코드의 일환인것으로 보임 : https://jakevdp.github.io/PythonDataScienceHandbook/04.13-geographic-data-with-basemap.html
  • 여기서는 conda 환경에서 basemap을 설치하라고 하는데 왠만하면 conda 환경 구축을 피하고 싶었다.
• numpy에서 투명도 조절하기 : black to transparent alpha image with numpy
  • https://www.kaggle.com/code/marcs1/black-to-transparent-alpha-image-with-numpy

2022년도 회고를 보아하니까 그 때랑 지금의 나랑 또 다르다는 생각이 들어서 신기했다.
그리고 작년 회고는 2021년도와 2022년도를 같이 묶어서 했기 때문에
시간이 이렇게 흘렀고 내 움직임은 대략적으로 이렇더라~ 라서
이번에는 키워드를 바탕으로 적어보겠다.

2023년도의 키워드

탈출과 이사, 처음 가는 길, 시각화, 커뮤니티

탈출과 이사

원래 이 키워드에는 '이사'만 있었다. 그런데 '탈출'을 아니 적을 수 없겠더라.
먼저 조직 개편과 팀 이동 2번의 이사가 있었다.
조직 개편을 하면서 탈출을 했고 팀 이동을 통해 이사를 했다.

탈출

탈출에 대해서 말하기 조심스럽지만 키워드로 넣은 이유가 있다.
2022년도에서 '좋은 엔지니어가 되기 위해서 고려해야 하는 점'을 적어놨던데
실은 2021년-2022년도에는 시간과 성과와 자원의 사용에 대해 많이 쫓겼다.

그런데! 조직 개편으로 그 분위기에서 벗어날 수 있었다는 말로 탈출에 대한 설명을 대신하겠다.
오죽하면 2023년도 초에는 다시 회사를 들어온 것 같은 느낌을 받았다.
아아니? 기술 개발 '연구'를 위해 이만큼 자원을 쓸 수 있다고???
다른 사람들이 보면 별거 아닌 것 처럼 보이겠지만 나에게는 감동이었다.

이사

그리고 직군도 좀 바뀌었다.
원래는 GIS엔지니어면서 데이터 엔지니어가 하는 일을 한번씩은 다 해봤는데(주의: 다 해봤다 해서 할줄 아는게 아니다)
지금은 데이터 전처리와 시각화, 분석에 초점을 맞추면서 연구개발쪽으로 많이 넘어왔다.
그래서 팀도 개발팀에서 연구팀으로 한번 더 바뀌게 되었다.
지금은 내 자신이 데이터 리서치 분야에 좀 가깝지 않나 생각이 들고 있는데, 내년은 또 어떻게 될 지 모르지.

처음 가는 길

완전히 조직도 팀도 바뀌니까 모든 것이 새로웠다. 여름쯤에 팀을 바꾸기 전에는 내가 여기 계속 있어도 되는가 하는 정체성의 혼란을 느꼈다.
으레 늘 하는 얘기지만 업무가 업무다 보니 내가 회사를 왕따시키는게 아닐까 하는 우스갯소리가 생각날 정도로 혼자다! 그리고 외로웠다!
2023년도부터는 회사에 천사같은 GIS 엔지니어분이 한분 더 합류하셔서 덜 고독하다.

그리고 또다른 고독은 문제 해결에서 찾아왔다.
내가 하려는 일에 대해 레퍼런스를 찾고 싶어도 찾을 수 없어서 머리를 정말 열심히 굴렸다.
예를 들어 A->E까지 가는 법을 알고 싶다 치자. 나는 A->B, B->C, C->D, D->E 혹은 그 이상으로 문제를 쪼갤 수 있을 만큼 쪼개고
해결법을 엮어 나가고 싶었는데 A->B가는 일도 보통이 아닐 때가 있다.
아, 이런 일을 하면서 문제를 잘 정의하는 법도 배웠지만 나는 매일이 불안했다.

이건 시간이 흐르고 팀의 소속감을 느끼면서 점차 잦아들었다.
아예 없어졌다고 할 수 없다. 나 혼자 결정을 내려야 한다면 나는 내 모든 결정에 대해서 판단하고 책임져야한다.
그런 쪽의 경각심은 늘 가지며 살고 있다.

시각화

처음 가는 길에서 나왔던 이야기는 사실 시각화 이야기다.
아니? 위성영상 기반 모델의 인퍼런스 결과를 지도 위에 시각화하라구요? 아아니 이걸 보고서로 만들라구요??!
시각화에 그치지 않고 분석으로 이어지는 일을 하다 보니 선례 삼을 만한게 없어서 눈뜬 장님이 된 기분이었다.
쓰다 보니까 내가 하반기 후반에 접어들기까지 정말 불안했나보다.

하반기에 내가 시각화한 결과물들이 전시 자료로 쓰이고 홍보 책자에 실리는 것을 보고 놀랐다.
아니 이렇게까지 반응이 좋다고!?
(적고 보니 아니 근데 진짜를 빼놓고 얘기 할 수 없는거니...)
거기에 더해 내가 작년에 한 일을 사람들이 기억하고 날 찾아줘서 놀랐다.
보안 때문에 뭐 한다고 말도 못하고 그렇게 끝난 프로젝트였기 때문이었다.

내가 만든 결과물이 영향력이 있구나를 확인할 수 있는 소중한 기회들이었다.

커뮤니티

원래 커뮤니티 활동에 소극적이었다.
ENFP고 사람은 좋아하는데 내 직무가 워낙 독특하니까 사람들 간의 연결고리를 찾기 힘들었다.
올해 여름에 페북 친구인 강성욱님이 나에게 먼저 손을 내밀어주셨다.
그렇다, 나는 누군가 나를 이끌면 가속도가 붙더라고.
올해가 다 가는 시점에서는 마침내 운영진에 합류할 수 있었다!
운영진 추천도 감사하게도 강성욱님이 나를 적극 추천해주셨다.
(링크 : https://k-devcon.com/staff)
커뮤니티에 가입한지 어연 반년이 지났는데, 이번에는 아이템으로 느낀 점들을 써보고자 한다.

커뮤니티 활동을 하면서 느낀 건 다음과 같다.

• 다른 사람들이 얼마나 치열하게 사는지
• 주말에 할거 없을 때 나와서 공부 좀 합시다
• 내가 겪는 고민을 똑같이 겪는 사람들이 많구나. 내가 괜한 걱정을 만든 게 아니었구나
• 2022년에 다짐했던 것 : 지식을 나눈다 -> 어떤 것이든 내가 나누는 것들이 지식이 되고 도움을 받는 사람들이 생긴다
• 내가 이렇게 사람들이랑 같이 어울리다 보면 그들에게도 조금씩 배울 것들이 생기고, 나도 그들에게 줄 수 있는게 생기는구나

운영진으로서 다음과 같은 생각을 하고 있다

• 이렇게 꾸준히 활동 하는(격주 토요일이다) 그룹에서 운영진으로 활동하는 것은 처음이라서 떨리고 긴장된다
• 커뮤니티 사람들이 여기서 편안함을 느끼고 조금씩 십시일반 해서 뭐라도 해봤으면 좋겠다. 당장 나도 그런 생각을 가지고 들어왔으니까
• 내년에는 조금 더 적극적인(소극적이여도 괜찮다. 소규모로 플래시몹 스터디를 짧게 짧게 할까 생각도 들고...) 활동을 해보자
• 리더쉽, 운영과 기획, 스트레스 관리에 대해서 배우자

내년 포부

2024년의 팀 내 로드맵도 내 로드맵도 대략 다 짜여있다.
이번 해도 공부할 것이 많고(분산 처리를 배워야 한다), 맡는 것도 의외로 많고(운영진과 운영진과 노사협의회), 2023년도에 했던 일을 좀 더 발전시켜야 한다. 너에게 확신을 가졌으면 좋겠다 혜미야. 불안했던 날들 방황했던 날들이 이제 완전히 끝이 난 건 아니지만 이제는 너를 믿어주는 사람들이 많으니까 네 자신부터 나를 많이 아꼈으면 좋겠다. 낙수가 주춧돌을 뚫듯 그렇게 잔잔하고 꾸준하며 우직한 힘을 내길 바라.
내년에도 잘 부탁해!

개요

해당 문서의 대상 독자(이하 서술 조건의 대부분에 속하신다면 대상이십니다)

• 급하게 ArcGIS로 데이터를 처리해야 함
• javascript보다는 python이 편함
• ArcGIS에 데이터는 올렸는데 편집이 필요하다. 근데 나 맥 써서...ArcGIS Pro로 대응을 할 수 없다. <- 글쓴이는 이 경우였음
• ArcGIS Online이 편하다.

해당 문서의 주요 키워드

• ArcGIS
• ArcGIS Online
• ArcGIS Python API

내용

1. ArcGIS API package를 통한 local 사용

전제 : ArcGIS Pro 또는 ArcGIS Server가 깔려 있는 컴퓨터

• 2021년까지만 해도 ESRI는 ArcGIS 제품이 해당 머신에 깔려 있는 경우에만 API를 사용하도록 허용했다.
• 2021년 부터는 새로운 방법을 지원한다고 하는데, 솔직히 뭐가 바뀌었는지는 모르겠다.(문서 작성 당시인 2022년 10월에 시도했지만 실패)
• 해당 방법은 기존 API 사용시 가장 기본으로 생각하는 경우이다.
• 두 제품군 모두 윈도우에서만 지원한다.
• arcgis (python) package가 미리 깔려있음을 전제함
• 다음의 두가지 방법이 있다(참고 : 문서)
  • python pakcage manager를 사용하라고 안내한다.
  • arcgis가 설치된 폴더 내에 위치한 Python Command Prompt를 켜서 conda명령어를 통해 python api를 깐다.

전제 : ArcGIS 제품군이 없는/설치 불가인 머신

(문서 작성 시 참 많은 시도를 해봤는데, 잘 되지 않았다.)

• conda를 사용해서 리눅스에도 깔 수 있다고 문서에 나와있지만 실패.
• install offline
  • 파일 다운로드 받아서 압축 풀기 > conda에서 offline 설치 가능하게 환경 설정 > 명령어
  • 공식 문서에는 The conda
    utility will pull all the arcgis package dependencies from the installed set of Anaconda libraries instead of searching the internet. 라고 되어있지만 cli로 테스트해보니까 제대로 import 되지 않아 포기
• docker install
  • 사실상 윈도우만 지원 (참고 : 문서)

2. ArcGIS Online을 통한 API 사용

전제 : ArcGIS Online 내 notebook 서버 사용

• arcgis package만 지원하는 서버 : 무료

  • 우리가 하는 일을 하려면 arcgis 서버에 접근해야 하는데, 접근 권한을 얻기 위해서는 arcgis package에서 지원하는 GIS class 를 통해 접근 권한을 얻어야 한다.

  • 예시

    from arcgis.gis import GIS
    gis = GIS() # anonymous connection to www.arcgis.com
    # Search for 'USA major cities' feature layer collection
    search_results = gis.content.search('title: USA Major Cities',
                                      'Feature Layer')

• 다른 옵션

  • arcpy 지원하는 서버 : 크레딧
  • arcpy + gpu resouces : 크레딧

ArcGIS Online을 통한 Python API 사용 예시

1. feature layer access

주의 : 아래 나오는 코드들은 모두 예시 코드로, 실제 접속 정보는 제외함. 데이터도 이하동문.

password = 'blablabla'
# jungham.map.arcgis.com은 해당 계정이 속한 organization의 주소
# source로 해당 계정이 접속할 수 있는 데이터에 접근 가능(contents)
source = GIS("https://jungham.maps.arcgis.com", f"{arcgis_user_account}", f"{password}", verify_cert=False)
# 이 레이어의 이름은 ham, 자원 유형은 'Feature Layer'이다 
sample_data_layer = source.content.search('ham',item_type='Feature Layer')

2. feature layer update

1) geometry 형태

{"x": -10617585.164750207,"y": 5160968.961654072,"spatialReference": {"wkid": 102100, "latestWkid": 3856}
}

2) feature 형태

아래는 layer를 이루는 feature 하나의 정보이다.

    {"geometry": 
    {"x": 13948339.367548944, "y": 4677527.208890018, "spatialReference": {"wkid": 4326}}, 
    "attributes": {"FID": 6, "type": 'candyshop', "content": "old-fashioned candyshop", "create_datetime": "2019-03-11 00:00:00"}
    }

3) 업데이트 예시

ArcGIS 내에서 데이터 구조를 이해하면 아래의 코드가 쉽게 이해된다.
OGC 표준에서는 보통 featurecollection > features > feature 이런 단계로 정의가 되는데
ArcGIS에서는 featurelayer > layer > featureset > feature 이런 구조로 데이터를 관리하고 있다.
그리고 Pandas의 DataFrame을 확장한 spatially enabled DataFrame(sedf)이라는 독특한 데이터 타입을 정의해서 featureset에 이 타입을 대응해서 사용하고 있다.

# spatially enabled Dataframe -> featureset -> update
# features in layer are treated as 'featureset'
sedf = pd.DataFrame.spatial.from_geodataframe(sample_event_df)
# adds, updates, deletes 등의 키워드가 있다. 
sample_data_layer[1].layers[0].edit_features(adds=sedf.spatial.to_featureset())

ArcGIS 내에서 데이터를 관리하는 구조와 그 타입

해당 내용은 Python API 사용 방법을 이해하는데 도움이 되었기에 사용 예제와 함께 적는다.

데이터 구조

• FeatureLayer > layer > featureset > feature
• 실제 서버에서 publish된 layer(feature layer를 예시로 들면) 아래에는 여러 layer가 포함되어 있을 수 있다. 그리고 하나의 layer는 featureset 내에 feature 데이터를 저장한다. featureset은 feature의 집합이다. feature들을 list로 가지고 있다.

시간 활성화

• 시간 활성화 메뉴얼 1 : https://support.esri.com/en/technical-article/000024836
• 시간 활성화 메뉴얼 2 : https://doc.arcgis.com/en/arcgis-online/reference/configure-time.htm

시간 활성화 메뉴얼에 대한 추가 설명

시간 활성화 메뉴는 메뉴얼을 봐도 찾기 힘들었다.
published feature layer 아래에 layer의 집합이 있다. 그러니까 feature layer는 layer의 집합이다.
예를 들면, 'the location of the opened shops' 라는 이름의 feaure layer가 있으면 그 아래에 'the location of the opened shops_0'이라는 point layer가 존재한다. 실제로 사용한 feature layer의 경우 {feature layer name}_0라는 이름을 가진 기하 레이어가 그 아래에 존재한다.

피처 레이어는 건물, 필지, 도시, 도로, 지진 진앙 등 유사한 지리적 피처를 그룹화한 것입니다. 피처는 포인트, 라인 또는 폴리곤(영역)일 수 있습니다.

문구 출처 : 피처 레이어

아래 첫번째 그림은 feature layer아래에 위치한 layer에 들어가서 시간 옵션을 활성화하는 예시이다.

그림 출처 : How To: Enable time on a layer in ArcGIS Online and create a web app with time animation

후기와 회고

• ArcGIS Python API를 사용하는데, 문서간의 연속성이 없어서 문서 간 이동을 통해 얻어낸 내용을 정리하고 싶었다.
• 좀 더 직설적으로 '구려!'라고 말하고 싶지만, 누구에게는 ArcGIS가 많은 일들을 해주고 있을 것이고, 실제로도 그렇기 때문에 그렇게는 적고 싶지 않았다. 문제는 존재하니까 남은건 푸는 것일 뿐. 다만 아래 한 마디만 남기겠다.
• 한국 서비스 지원 센터에 문의를 했으나...
• ArcGIS 한국인 사용자 중 python api를 어떻게 써야 할지 감이 잡히지 않는다면, 데이터를 불러오는 간단한 예제부터 따라해보라고 말하고 싶다. 이 글이 그런 분들에게 도움이 되길 바라는 마음에서 작성했다.
• 정리하다 보니, feature layer crud 예제까지 작성해도 될 것 같았다. 그것은...! 시간이 되면 이 글에다가 추가하겠다.

개요

• qgis 3.22버전에서 새롭게 annotation layer라는 layer가 새로 생겼다.
  qgis에서 annotation layer가 포함된 qgs를 열면 다음과 같이 보인다. raster와 qgs가 한 쌍으로 움직인다.
  

  

• 이 layer는 벡터 데이터가 아니며, qgis project file(qgs)에 xml로 도형을 기록한다
• 어떻게 annotation의 geometry를 추출할 수 있는가?

참고 자료

• qgis annotation layer란?
• 3.22에는 annotation layer을 대응할 수 있는 tool bar가 있는 것으로 보이나, 내가 쓰는 QGIS 3.26 버전에서는 찾을 수 없었다.
  • (수정 : 22.11.18) macos QGIS 3.28 버전 기준, View > Toolbar > annotation toolbar(주석 툴바)로 추가하시면 다음의 툴바를 확인할 수 있습니다. 글 작성 시점으로는 3.29까지 나왔습니다. 

• annotation layer 정보는 qgs 내에 xml로 기술되어있다.
  

  

  
  위 그림에서 볼 수 있듯, qgs내에 기술된 annotation layer와 qgs가 바라보는 raster layer 두 레이어가 뜨나, qgz를 압축해제해서 얻은 결과는 qgs와 링크된 raster 데이터 둘 뿐이었다.

코드

import shapely.wkt
import xml.etree.ElementTree as elemTree
from osgeo import ogr
import geopandas as gpd

'''
box, polygon, mask
'''
'''
parse_annotation은 여기서 qgs 내에 xml으로 기술되어있는 프로젝트 정보 중 'layer'단위의 xml을 파싱한다. 
return : annotation 파일 내의 geometry를 모아서 shapely.geometry 타입으로 list를 반환
'''
def parse_annotation(elem) : 
    items = elem.findall('.//items/item[@wkt]')
    wkts = [x.get('wkt') for x in items]
    geoms = list()
    for wkt in wkts : 
        if 'Curve' in wkt :
            try : 
                geom = ogr.CreateGeometryFromWkt(wkt)
                geom_approx = geom.GetLinearGeometry()
                geom_shapely = shapely.wkt.loads(geom_approx.ExportToWkt())
                geoms.append(geom_shapely)
            except : 
                print(wkt)
        else : 
            geom_shapely = shapely.wkt.loads(wkt)
            geoms.append(geom_shapely)
    return geoms

'''
path : 파일의 위치
return : annotation 파일 내의 전체 geometry를 list로 반환 
'''

def parse_qgs_for_annotation(path) : 
    tree = elemTree.parse(path)
    root = tree.getroot()

    main_annotation_layer = root.findall('.//main-annotation-layer')[0]
    extra_annotation_layers = root.findall(".//projectlayers/*[@type='annotation']")

    geoms = list()
    geoms.extend(parse_annotation(main_annotation_layer))
    for extra_annotation_layer in extra_annotation_layers :
        geoms.extend(parse_annotation(extra_annotation_layer))
    return geoms

주의사항

CurvePolygon 때문에 ogr 모듈을 썼다. esri shapefile format에서는 geometry에서 curve 타입을 지원하지 않는다고 했기 때문에, shapefile 말고 geojson으로 내보내길

용례

tree = elemTree.parse(path)
root = tree.getroot()
srid = root.find('./projectCrs//srid')
srid = int(srid.text)

res = parse_qgs_for_annotation(path)
gdf = gpd.GeoDataFrame(geometry=res)
gdf = gdf.set_crs(f'epsg:{srid}')

gdf.to_file("test.geojson", driver='GeoJSON')