쉽게 시작하는 pygis

쉽게 시작하는 PyGIS, 총 6개의 레슨으로 이루어져있는 수업자료를 기반으로 진행합니다.

원본 문서 : https://automating-gis-processes.github.io/CSC18/

이번 문서는 Lesson1을 다룹니다. 실습 문서 링크

Intro

설치

PyGIS를 시작하기 위해 설치해야 패키지가 있습니다. 오늘 실습에서 사용하는 패키지는 아래 리스트면 됩니다.

import shapely

PyGIS로 다음과 같은 것을 할 수 있습니다.

1. 다양한 파일 형식에서 / 다른 공간 형식으로 공간 데이터 읽기 / 쓰기
2. 서로 다른 투영기법에 대해 다루기
3. 서로 다른 기하 연산과 공간 질의 수행
4. 주소를 위경도로 변환 (혹은 그 반대) 즉, 지오 코딩
5. 다른 기준에 따라 데이터를 재분류
6. Python으로 OpenStreetMap에서 데이터를 쉽게 가져 오는 방법 알아보기
7. 파이썬에서 래스터 처리의 기초 배우기
8. 지리 데이터를 시각화 및 (대화형)맵을 표현

(사실상 만능이네요!) 무엇보다 제일 좋은 이유는 무료이기 때문이죠. 그리고 파이썬이기 때문!

그리고 다양한 GIS툴들이 있지만 해당 내용의 자세한 설명은 아래의 문서를 참조하시면 됩니다!

Data analysis & visualization:

• Numpy –> Fundamental package for scientific computing with Python
• Pandas –> High-performance, easy-to-use data structures and data analysis tools
• Scipy –> A collection of numerical algorithms and domain-specific toolboxes, including signal processing, optimization and statistics
• Statsmodels –> Statistical models for Python
• Scikit-learn –> Machine learning for Python (classification, regression, clustering, etc.)
• Matplotlib –> Basic plotting library for Python
• Seaborn –> Statistical data visualization
• Bokeh –> Interactive visualizations for the web (also maps)
• Plotly –> Interactive visualizations (also maps) for the web (commercial - free for educational purposes)
• Dash –> Building analytical web applications with Python (no Javascript required)and GIS Software Library for Python.

GIS:

• GDAL –> Fundamental package for processing vector and raster data formats (many modules below depend on this). Used for raster processing.
• Geopandas –> Working with geospatial data in Python made easier, combines the capabilities of pandas and shapely.
• Shapely –> Python package for manipulation and analysis of planar geometric objects (based on widely deployed GEOS).
• Fiona –> Reading and writing spatial data (alternative for geopandas).
• Pyproj –> Performs cartographic transformations and geodetic computations (based on PROJ.4).
• PyCRS –> Working eaily with different CRS specifications (EPSG, ESRI, Proj4)
• Pysal –> Library of spatial analysis functions written in Python.
• Geopy –> Geocoding library: coordinates to address <-> address to coordinates.
• GeoViews –> Interactive Maps for the web.
• Geoplot –> High-level geospatial data visualization library for Python.
• GeoNotebook –> Desktop GIS-like environment for visualizing and interacting with spatial data using Python (based on Jupyter Notebooks)
• OSMnx –> Python for street networks. Retrieve, construct, analyze, and visualize street networks from OpenStreetMap
• Networkx –> Network analysis and routing in Python (e.g. Dijkstra and A* -algorithms), see this post.
• Cartopy –> Make drawing maps for data analysis and visualisation as easy as possible.
• Scipy.spatial –> Spatial algorithms and data structures.
• Rtree –> Spatial indexing for Python for quick spatial lookups.
• Rasterio –> Clean and fast and geospatial raster I/O for Python.
• Rasterstats –> A module for summarizing geospatial raster datasets based on vector geometries (e.g. conduct zonal statistics).
• RSGISLib –> Remote Sensing and GIS Software Library for Python.

Lesson1

이번 레슨에서는 공간 데이터 모델 (기하학적 객체)과 순수 Python에서 공간 데이터를 조작하는 방법에 대한 몇 가지 기본 사항을 소개합니다.

1. 기하학적 객체 - 공간 데이터 모델
   • 기하학적 객체 및 모양새 모듈 개요
   • 포인트 객체
   • 선형 문자열 객체
   • 다각형 객체
   • 지오메트리 컬렉션(선택 사항) 2.실습 1 : 기하학적 객체 작업 3.실습 1 정답

학습 목표

1. 파이썬에서 GIS를 수행하는 데 사용할 수있는 도구에 대한 아이디어가 있습니다.
2. conda를 사용하여 Python 패키지를 설치하는 방법을 알아야합니다.
3. 어떤 종류의 기하학적 객체가 사용 가능한지 알고있다.
4. Shapely를 사용하여 다른 종류의 지오메트리를 작성하는 방법을 알아야합니다.
5. 파일에서 좌표를 읽고 이를 바탕으로 점을 만드는 법을 알아야합니다.

기하학적 객체 및 모양새 모듈 개요

shaply 모듈은 점,선,다각형의 기하하적 객체를 만들고 작업하는 패키지입니다.

기하 객체는 좌표 튜플로 구성됩니다.

Point-

object는 공간상의 단일 점을 나타냅니다. 점은 2 차원 (x, y) 또는 3 차원 (x, y, z) 일 수 있습니다.

LineString

-object (즉, 선)는 선을 형성하기 위해 서로 조인 된 점들의 순서를 나타냅니다. 따라서 선은 적어도 두 개의 좌표 튜플 목록으로 구성됩니다.

Polygon

-object는 외부 링을 형성하는 적어도 세 개의 좌표 튜플 목록과 구멍 다각형 목록 (가능한)으로 구성된 채워진 영역을 나타냅니다.

기하학적 객체의 집합 (예 : 다각형이있는 다각형)을 가질 수도 있습니다.

MultiPoint

-object는 포인트 모음을 나타내며 좌표 튜플 목록으로 구성됩니다.

MultiLineString

-object는 라인 모음을 나타내며 라인 계열 시퀀스 목록으로 구성됩니다.

MultiPolygon

-object는 외부 링과 (가능한) 홀리스트 튜플을 구성하는 다각형과 같은 시퀀스의 목록으로 구성된 폴리곤 콜렉션을 나타냅니다.

Point

코드로 표현해보면 다음과 같습니다.

  # 1. shapely 패키지 호출
  from shapely.geometry import Point, LineString, Polygon

  # 2. 점Point 지리 객체 생성하기
  point1 = Point(2.2, 4.2)
  point2 = Point(7.2, -25.1)
  point3 = Point(9.26, -2.456)
  point3D = Point(9.26, -2.456, 0.57)

  # 3. 점의 타입 확인
  point_type = type(point1)

  # 4. 결과값
  In [7]: print(point1)
  POINT (2.2 4.2)

  In [8]: print(point3D)
  POINT Z (9.26 -2.456 0.57)

  In [9]: print(type(point1))
  <class 'shapely.geometry.point.Point'>

  # 5. 좌표값 추출
  In [10]: point_coords = point1.coords

  # 6. 좌표값 타입
  In [11]: type(point_coords)
  Out[11]: shapely.coords.CoordinateSequence

  # 7. 좌표값에서 x,y 속성
  In [12]: xy = point_coords.xy

  # 8. x좌표값
  In [13]: x = point1.x

  # 9. y좌표값
  In [14]: y = point1.y

  # 10. 좌표결과
  In [15]: print(xy)
  (array('d', [2.2]), array('d', [4.2]))

  In [16]: print(x)
  2.2

  In [17]: print(y)
  4.2

앞으로도 자주 쓰일 두 점 사이의 거리를 구하는 방법은 무척 간단합니다. 점1.distance(점2) 끝!

# Calculate the distance between point1 and point2
In [18]: point_dist = point1.distance(point2)

In [19]: print("Distance between the points is {0:.2f} decimal degrees".format(point_dist))
Distance between the points is 29.72 decimal degrees

LineString

LineString -object를 생성하는 것은 Point가 생성되는 방법과 매우 유사합니다. 이제는 하나의 좌표 튜플을 사용하는 대신 매끈한 Point 객체 목록이나 좌표 튜플을 사용하여 선을 구성 할 수 있습니다.

# 1. 점을 이은 선을 만들기
In [20]: line = LineString([point1, point2, point3])

# 2. 튜플형태의 좌표값을 직접 입력
In [21]: line2 = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])

# 3. 결과값(동일)
In [22]: print(line)
LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)

In [23]: print(line2)
LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)

In [24]: type(line)
Out[24]: shapely.geometry.linestring.LineString

LineString 속성 및 함수

LineString -object는 많은 유용한 내장 속성과 기능을 가지고 있습니다. 예를 들어 LineString (선)의 좌표 또는 길이를 추출하고, 선의 중심을 계산하고, 특정 거리에서 선을 따라 점을 작성하고, 선에서 지정된 점까지의 가장 가까운 거리를 계산하고 형상을 단순화하는 것이 가능합니다. Shapely 문서 의 전체 기능 목록을 참조하시면 됩니다.

# 1. 라인에서 xy좌표 추출하기
In [25]: lxy = line.xy

In [26]: print(lxy)
(array('d', [2.2, 7.2, 9.26]), array('d', [4.2, -25.1, -2.456]))

numpy array형태로 추출된것을 확인할 수 있습니다. 첫번째 array는 X좌표들을, 두번째 array Y좌표들을 가진것을 각각 독립적으로 가집니다. 그렇기 때문에 아래처럼 특정 좌표만을 뽑을 수 있습니다.

# 1. X좌표 추출
In [27]: line_x = lxy[0]

# 2. index로 y좌표 추출
In [28]: line_y = line.xy[1]

In [29]: print(line_x)
array('d', [2.2, 7.2, 9.26])

In [30]: print(line_y)
array('d', [4.2, -25.1, -2.456])

또, 라인의 길이. 라인의 중심(centroid) 정보 또한 가져올 수 있습니다.

# 1. 라인의 길이
In [31]: l_length = line.length

# 2. 라인의 중심점
In [32]: l_centroid = line.centroid

# 3. 중심점의 타입
In [33]: centroid_type = type(l_centroid)

# 4. 결과 출력
In [34]: print("Length of our line: {0:.2f}".format(l_length))
Length of our line: 52.46

In [35]: print("Centroid of our line: ", l_centroid)
Centroid of our line:  POINT (6.229961354035622 -11.89241115757239)

In [36]: print("Type of the centroid:", centroid_type)
Type of the centroid: <class 'shapely.geometry.point.Point'>

Polygon

다각형 객체는 Point 및 LineString이 작성된 방법과 동일한 논리를 계속하지만 Polygon 객체는 입력으로 좌표 튜플만 허용합니다. 다각형에는 최소한 세 개의 좌표 튜플이 필요합니다.

# 1. 세 개의 좌표로 다각형 만들기
In [37]: poly = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])

# 2. 점 객체로 생성하기(x,y 좌표계가 반드시 있어야 함)
In [38]: poly2 = Polygon([[p.x, p.y] for p in [point1, point2, point3]])

# 3. 지리 타입을 문자열로 접근가능
In [39]: poly_type = poly.geom_type

# 4. 지리 타입을 다른 포맷으로 저장
In [40]: poly_type2 = type(poly)

# 4. 결과값
In [41]: print(poly)
POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))

In [42]: print(poly2)
POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))

In [43]: print("Geometry type as text:", poly_type)
Geometry type as text: Polygon

In [44]: print("Geometry how Python shows it:", poly_type2)
Geometry how Python shows it: <class 'shapely.geometry.polygon.Polygon'>

Polygon에는 좌표 주위에 이중 괄호가 있습니다. Polygon에는 내부에 구멍이있을 수도 있기 때문입니다. Polygon -object의 도움으로, 외부 좌표와 내부 좌표 (선택 사항)를 사용하여 Polygon을 구성 할 수 있습니다. 내부 좌표가 Polygon 내부에 구멍을 만듭니다.

Help on Polygon in module shapely.geometry.polygon object:
class Polygon(shapely.geometry.base.BaseGeometry)
 |  A two-dimensional figure bounded by a linear ring
 |
 |  A polygon has a non-zero area. It may have one or more negative-space
 |  "holes" which are also bounded by linear rings. If any rings cross each
 |  other, the feature is invalid and operations on it may fail.
 |
 |  Attributes
 |  ----------
 |  exterior : LinearRing
 |      The ring which bounds the positive space of the polygon.
 |  interiors : sequence
 |      A sequence of rings which bound all existing holes.
Let’s create a Polygon with a hole inside

다각형 내부를 그려봅시다.

# 1. 먼저 외부를 정의합니다.
In [45]: world_exterior = [(-180, 90), (-180, -90), (180, -90), (180, 90)]

# 2. 내부 구멍을 정의 합니다.
주의 : 다중 구멍을 만들때는 리스트 형태로 넣어줘야 합니다.
In [46]: hole = [[(-170, 80), (-170, -80), (170, -80), (170, 80)]]

# 3. 구멍이 없는 공간을 그립니다.
In [47]: world = Polygon(shell=world_exterior)

# 4. 이제 구멍이 있는 공간을 그려봅시다.
In [48]: world_has_a_hole = Polygon(shell=world_exterior, holes=hole)

# 5. 결과값
In [49]: print(world)
POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90))

In [50]: print(world_has_a_hole)
POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90), (-170 80, -170 -80, 170 -80, 170 80, -170 80))

In [51]: type(world_has_a_hole)
Out[51]: shapely.geometry.polygon.Polygon

다각형 속성 및 함수

Polygon의 면적, 중심, 경계 상자, 외부 및 외부 길이와 같이 실제로 유용한 다른 속성에 다시 액세스 할 수 있습니다.

# 1. 다각형의 중심
In [52]: world_centroid = world.centroid

# 2. 다각형의 면적
In [53]: world_area = world.area

# 3. 다각형의 둘레 (i.e. bounding box)
In [54]: world_bbox = world.bounds

# 4. 내부가 빈 다각형 외부 객체
In [55]: world_ext = world.exterior

# 5. 내부가 빈 다각형 외부 길이
In [56]: world_ext_length = world_ext.length

결과를 봅시다.

In [57]: print("Poly centroid: ", world_centroid)
Poly centroid:  POINT (-0 -0)

In [58]: print("Poly Area: ", world_area)
Poly Area:  64800.0

In [59]: print("Poly Bounding Box: ", world_bbox)
Poly Bounding Box:  (-180.0, -90.0, 180.0, 90.0)

In [60]: print("Poly Exterior: ", world_ext)
Poly Exterior:  LINEARRING (-180 90, -180 -90, 180 -90, 180 90, -180 90)

In [61]: print("Poly Exterior Length: ", world_ext_length)
Poly Exterior Length:  1080.0

지오메트리 컬렉션 (선택 사항)

이 부분은 필수는 아니지만 지오메트리 모음 및 일부 특수 기하학적 객체 (예 :경계 상자)의 구성 및 사용과 관련된 유용한 정보가 포함되어 있습니다.

어떤 경우에는 단일 피처 아래에 여러 선이나 다각형을 저장하는 것이 좋습니다 (예 : Shape 파일의 단일 행이 두 개 이상의 선 또는 다각형 객체를 나타냄). MultiPoint 객체를 사용하여 곡선 컬렉션을 구현하고 MultiLineString 객체를 사용하여 곡선 컬렉션을 만들고 MultiPolygon 객체로 표면 컬렉션을 구현합니다. 이러한 컬렉션은 계산 상 중요하지 않지만 특정 종류의 기능을 모델링하는 데 유용합니다. Y 모양의 선 형상 (예 : 도로) 또는 여러 다각형 (예 : 같은 섬)은 MultiLineString 또는 MultiPolygon을 사용하여 전체적으로 멋지게 표현할 수 있습니다. 최소 경계 상자 만들기 및 시각화 예를 들어 데이터 포인트 주변은 많은 목적 (예 : 데이터의 범위를 이해하려고 시도)에서 실제로 유용한 기능입니다. 여기서 Shapely를 사용하여 데이터 포인트를 만드는 방법을 보여줍니다.

In [62]: from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon, box

# 다중점을 생성하기.
In [63]: multi_point = MultiPoint([point1, point2, point3])

# 다중점은 튜플로도 가능합니다.
In [64]: multi_point2 = MultiPoint([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])

# 두개의 선으로 이루어진 다중선 만들기
In [65]: line1 = LineString([point1, point2])
In [66]: line2 = LineString([point2, point3])
In [67]: multi_line = MultiLineString([line1, line2])

# 동일하게 다중다각형도 생성합니다. 이때 동서로 나눕니다.
# 서쪽 지역
In [68]: west_exterior = [(-180, 90), (-180, -90), (0, -90), (0, 90)]

# 서쪽 지역의 구멍부분을 만듭니다.
In [69]: west_hole = [[(-170, 80), (-170, -80), (-10, -80), (-10, 80)]]

# 다각형을 생성
In [70]: west_poly = Polygon(shell=west_exterior, holes=west_hole)

#이제, 동쪽 지역의 경계박스를 하나 만듭니다. 이떄 좌측하단과 우측상단의 경계 좌표를 정해야합니다.
In [71]: min_x, min_y = 0, -90
In [72]: max_x, max_y = 180, 90

# box() function을 활용해서 동쪽 다각형을 만듭니다.
In [73]: east_poly_box = box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)

# 다중다각형을 다음과 같이 만듭니다.
In [74]: multi_poly = MultiPolygon([west_poly, east_poly_box])

In [75]: print("MultiPoint:", multi_point)
MultiPoint: MULTIPOINT (2.2 4.2, 7.2 -25.1, 9.26 -2.456)

In [76]: print("MultiLine: ", multi_line)
MultiLine:  MULTILINESTRING ((2.2 4.2, 7.2 -25.1), (7.2 -25.1, 9.26 -2.456))

In [77]: print("Bounding box: ", east_poly_box)
Bounding box:  POLYGON ((180 -90, 180 90, 0 90, 0 -90, 180 -90))

In [78]: print("MultiPoly: ", multi_poly)
MultiPoly:  MULTIPOLYGON (((-180 90, -180 -90, 0 -90, 0 90, -180 90), (-170 80, -170 -80, -10 -80, -10 80, -170 80)), ((180 -90, 180 90, 0 90, 0 -90, 180 -90)))

지오메트리 컬렉션 속성

# 다중다각형의 볼록한 부분(Convex Hull)을 구할 수 있습니다.--> https://en.wikipedia.org/wiki/Convex_hull
In [79]: convex = multi_point.convex_hull

# 다중선에서 몇개의 라인이 들어가 있는지 구할 수 있습니다.
In [80]: lines_count = len(multi_line)

# 다중다각형의 면적을 구할 수 있습니다.
In [81]: multi_poly_area = multi_poly.area

# 인덱스를 통해 특정 다각형에 접근할 수 있습니다. 서쪽 지역을 추출 해봅니다.
In [82]: west_area = multi_poly[0].area

# 유효한 "MultiPolygon"이 있는지 확인할 수 있습니다. MultiPolygon은 개별 폴리곤이 서로 교차하지 않으면 유효한 것으로 간주됩니다. 여기서 폴리곤은 공통의 '0-자오선'을 가지기 때문에 유효한 폴리곤
이 없어야 합니다. 이것은 데이터에서 토폴로지 오류를 찾으려고 할 때 매우 유용한 정보가 될 수 있습니다 .
# 각각의 다각형이 유효하다면 유효한 다중다각형인지 확인 할 수 있습니다.
In [83]: valid = multi_poly.is_valid

결과를 확인 해보면

In [84]: print("Convex hull of the points: ", convex)
Convex hull of the points:  POLYGON ((7.2 -25.1, 2.2 4.2, 9.26 -2.456, 7.2 -25.1))

In [85]: print("Number of lines in MultiLineString:", lines_count)
Number of lines in MultiLineString: 2

In [86]: print("Area of our MultiPolygon:", multi_poly_area)
Area of our MultiPolygon: 39200.0

In [87]: print("Area of our Western Hemisphere polygon:", west_area)
Area of our Western Hemisphere polygon: 6800.0

In [88]: print("Is polygon valid?: ", valid)
Is polygon valid?:  False