반응형
반응형

1. QGIS 파이썬 콘솔로 작업하기

  • QGIS에서 파이썬 콘솔 실행(CTRL + ALT + P)
    • 아래의 텍스트 편집기가 표시 되며 파이썬 코드 실행 가능

2. 레이어 로딩하기

  • 파이썬 콘솔 창에 아래 코드 입력
  • 2가지 방법이 있음(QgsProject, iface)

2.1 방법1

  • QgsProject 의 addMapLayers method 사용
link = QgsVectorLayer(r'D:\data\[2021-04-05]NODELINKDATA\MOCT_LINK.shp', 'mLink','ogr')
node = QgsVectorLayer(r'D:\data\[2021-04-05]NODELINKDATA\MOCT_NODE.shp', 'mNode','ogr')
QgsProject.instance().addMapLayers([link, node])

  • QgsVectorLayer : path 파일을 통해 Vector Layer 객체를 생성
    • path : 파일 경로 지정(메모리일 경우 별도 포맷으로 입력)
    • 레이어명칭 : 레이어 대표 명칭 입력
    • provider : 동작을 수행할 provider 입력
  • addMapLayers : vector Layer객체를 Map layer에 추가한다.

2.2 방법2

  • iface 사용
ilink = iface.addVectorLayer(r'D:\data\[2021-04-05]NODELINKDATA\MOCT_LINK.shp', 'ifaceLink','ogr')
inode = iface.addVectorLayer(r'D:\data\[2021-04-05]NODELINKDATA\MOCT_NODE.shp', 'ifaceNode','ogr')

3. 레이어 조회, 삭제

3.1 레이어 조회

  • QgsProject 클래스의 mapLayers method 사용
    • EX) 4개의 레이어 객체가 출력됨
QgsProject.instance().mapLayers()
{'ifaceLink_32ac8313_af37_4c95_a181_274d20739ec1': <QgsVectorLayer: 'ifaceLink MOCT_LINK' (ogr)>,
 'ifaceNode_848b785e_3b87_4741_96c5_d02330f04282': <QgsVectorLayer: 'ifaceNode MOCT_NODE' (ogr)>,
 'mLink_4286530d_ef44_4b31_9326_af0d07f501c2'    : <QgsVectorLayer: 'mLink' (ogr)>,
 'mNode_aee8efc1_f39c_4955_807d_a953a0cb3e6d'    : <QgsVectorLayer: 'mNode' (ogr)>}

3.2 레이어 삭제

  • QgsProjet 클래스의 removeMapLayer method 사용
    • 인자로 상기 결과의 id값을 직접 전달
    • 또는 layer명.id() 를 인자로 전달
# 1. id값을 직접 전달
QgsProject.instance().removeMapLayer('ifaceLink_32ac8313_af37_4c95_a181_274d20739ec1')

# 2. layer명.id() 를 인자로 전달
QgsProject.instance().removeMapLayer(node.id()))

# 3. 다시 레이어 조회
QgsProject.instance().mapLayers()
{'ifaceNode_848b785e_3b87_4741_96c5_d02330f04282': <QgsVectorLayer: 'ifaceNode MOCT_NODE' (ogr)>,
 'mLink_4286530d_ef44_4b31_9326_af0d07f501c2'    : <QgsVectorLayer: 'mLink' (ogr)>}

4. 레이어 속성

4.1 레이어 좌표계 확인

  • 좌표계 대표 명칭 확인
crs = link.crs()
crs.description()
'ITRF2000_Central_Belt_60'
  • 좌표계 WKT로 표현
crs.toWkt()
'PROJCS[
    "ITRF2000_Central_Belt_60"
    ,GEOGCS["ITRF2000"
        ,DATUM["International_Terrestrial_Reference_Frame_2000"
            ,SPHEROID["GRS 1980",6378137,298.257222101]
            ,AUTHORITY["EPSG","6656"]
        ]
        ,PRIMEM["Greenwich",0]
        ,UNIT["Degree",0.0174532925199433]
    ]
    ,PROJECTION["Transverse_Mercator"]
    ,PARAMETER["latitude_of_origin",38]
    ,PARAMETER["central_meridian",127]
    ,PARAMETER["scale_factor",1]
    ,PARAMETER["false_easting",200000]
    ,PARAMETER["false_northing",600000]
    ,UNIT["metre",1,AUTHORITY["EPSG","9001"]]
    ,AXIS["Easting",EAST]
    ,AXIS["Northing",NORTH]
]'

4.2 레이어 상자 경계 얻기

  • string 형태로 레이어 상자 경계 얻기
extent = link.extent()
extent.toString()
'101766.5375999962707283,67516.4931999971158803 : 546275.7760999957099557,665746.5176000001374632'
  • Polygon(WKT) 형태로 레이어 상자 경계 얻기
extent.asWktPolygon()
'POLYGON((101766.5375999962707283 67516.49319999711588025, 546275.77609999570995569 67516.49319999711588025, 546275.77609999570995569 665746.51760000013746321, 101766.5375999962707283 665746.51760000013746321, 101766.5375999962707283 67516.49319999711588025))'
  • 상하좌우 경계 얻기
extent.xMinimum()
extent.xMaximum()
extent.yMinimum()
extent.yMaximum()
101766.53759999627
546275.7760999957
67516.49319999712
665746.5176000001

4.3 레이어 개수 얻기

  • featureCount
link.featureCount()
528232

4.4 레이어의 특정 feature 얻기

  • 첫 번째 feature의 geometry 조회
item = link.getFeatures() # feature iterator 생성
feat = QgsFeature() # 빈 feature 생성
item.nextfeatrue(feat) # feat에 첫 번째 feature 값 할당
feat.geometry() # 첫번 째 feature의 geometry 조회
<QgsGeometry: MultiLineString ((245889.2084229375468567 602540.10316239262465388, 245884.52438196362345479 602550.70766398275736719, 245880.58459971970296465 602562.44166173494886607, 245877.37697640058468096 602577.55608172016218305, 245874.80402102915104479 602590.92314239405095577, 245870.08166344027267769 602608.65557601337786764, 245866.1486026662751101 602619.1390588756185025, 245861.72407481138361618 602627.99418399820569903, 245850.4775304970680736 602649.69330804911442101, 245843.77653209323761985 602663.28828775370493531, 245840.21593625328387134 602674.27399252017494291, 245836.27413596049882472 602686.38314351008739322, 245833.96002093932474963 602698.12587899703066796, 245831.74675705164554529 602714.37113853613846004, 245830.48617608763743192 602739.50042814784683287, 245829.72781690946430899 602764.25725230504758656, 245829.78834916753112338 602799.52310311701148748, 245829.1886519008257892 602818.02802740619517863, 245827.46483578640618362 602836.27679859660565853, 245825.86203727853...>
  • vector type 확인
feat.geometry().type()
1

4.5 필드 정보 조회

  • 속성 정보 개수 조회
feat.fields().count()
17
  • 특정 속성정보 이름/Type 조회
feat.fields()[2].name()
feat.fields()[2].typeName()
'T_NODE'
'String'
  • 속성값 조회
    • 속성 이름, 속성 번호로 조회 가능
feat['T_NODE']
feat[2]
'2630076901'
'2630076901'

4.6 반복적으로 속성 조회

  • 도로명칭이 '화악산로' 인 도로 조회
    • 개수가 많아 5개 만 찾고 for문 종료
idx = 0
for f in link.getFeatures():
    if f["ROAD_NAME"] == '화악산로':
        print("OK idx = ", idx)
        idx = idx+1
    if(idx == 5):
        break
OK idx =  0
OK idx =  1
OK idx =  2
OK idx =  3
OK idx =  4

5. 신규 레이어 만들기(Point 예제)

  • 메모리에 레이어 생성
theLayer=QgsVectorLayer('Point?crs=epsg:4326','SomePoints','memory')
  • 속성 필드 추가
    • ID 와 NAME Field 추가
theFeatures = theLayer.dataProvider() # Feature를 다루는 class 할당
theFeatures.addAttributes([QgsField("ID", QVariant.Int),QgsField("Name", QVariant.String)])
  • 레이어에 feature 추가 하기
p=QgsFeature()   # 빈 feature 생성
point = QgsPointXY(-106.3463, 34.9685) # 임시 Geometry 클래스 생성
p.setGeometry(QgsGeometry.fromPointXY(point)) # 임시 Geometry를 p feature에 set
p.setAttributes([123,"Paul"]) # 속성 정보를 p feature에 set
theFeatures.addFeatures([p]) # Features class에 p를 add
theLayer.updateExtents() # 범위를 업데이트
theLayer.updateFields() # 속성정보를 업데이트
QgsProject.instance().addMapLayers([theLayer]) # 레이어 로딩

6. PostGIS와의 연동(Polygon 예제)

  • Sample 데이터를 postgres sql에 저장하고 실습 진행

    • Sample을 위해 행정 경계 데이터를 여기 에서 다운
    • 해당 데이터의 좌표계를 EPSG4326으로 변경 후 postgresql에 저장
  • postgresql 서버에 연결하여 polygon 데이터 조회하기


# 1. postgresql 서버에 연결하여 polygon 데이터 조회하기
import psycopg2
connection = psycopg2.connect(database='geospatial', user='postgres', password='qkrtkddus!1')
cursor = connection.cursor()
cursor.execute("select *, ST_AsTexT(geom) from tl_scco_ctprvn")
c = cursor.fetchall()

# 2. QGIS 레이어 및 Feature 생성(속성 필드 추가)
sigungu = QgsVectorLayer('Polygon','Sigungu','memory')
sigunguFeatures = sigungu.dataProvider()
sigunguFeatures.addAttributes([QgsField("ID", QVariant.Int), QgsField("Name",QVariant.String)])

# 3. QGIS 레이어에 postgresql 에서 조회해온 데이터 삽입하기
for poly in c:
  g=QgsGeometry.fromWkt(poly[5])  # geometry 저장
  p=QgsFeature()      # 빈 feature 생성
  p.setGeometry(g)    # feature 에 geometry setting
  p.setAttributes([poly[1], str(poly[3])]) # feature 에 속성 setting
  sigunguFeatures.addFeatures([p]) # features 변수에 feature 추가
  sigungu.updateExtents() # 범위를 업데이트
  sigungu.updateFields() # 속성정보를 업데이트
  QgsProject.instance().addMapLayers([sigungu]) # 레이어 로딩

7. Feature 추가, 편집, 삭제, 조회

  • 앞선 2. 레이어 로딩하기에서 로딩한 Node 데이터로 예제 진행
  • 레어어에 가능한 작업 리스트 확인
node.dataProvider().capabilitiesString()
'객체 추가, 객체 삭제, 속성 값 변경, 속성 추가, 속성 삭제, 속성 이름 바꾸기, 공간 인덱스 생성, 속성 인덱스 생성, ID로 빠른 객체 접근, 도형 변경'

7.1 Feature 추가

  • 새로운 Node 1개 추가
feat = QgsFeature(node.fields())
feat.setAttribute('NODE_ID', 99999999999)
feat.setAttribute('NODE_TYPE', 1)
# 또는 위의 코드를 아래와 같이 한줄로 표현 가능
# feat.setAttributes([999999999999, 1 , ...속성값들])
feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(221689.09,462751.59)))
node.dataProvider().addFeatures([feat])
추가 전 추가 후

7.2 Feature 삭제

  • 차선이 1개인 링크 삭제
deleteList = []
for x in link.getFeatures():
  if x["LANES"] == 1:
    deleteList.append(x.id())
link.dataProvider().deleteFeatures(deleteList)
삭제 전 삭제 후

7.3 Feature 편집

  • '서울외곽순환고속도로' 링크를 '수도권제1순환고속도로'로 명칭 변경
for x in link.getFeatures():
  if x['ROAD_NAME'] == '서울외곽순환고속도로':
    link.dataProvider().changeAttributeValues({x.id():{7:'수도권제1순환고속도로'}})
변경 전 변경 후

7.4 Feature 선택(수식 이용)

  • 도로번호가 1이면서 도로명이 '경부고속도로인 링크 선택하기
link.selectByExpression("ROAD_NO=1 and ROAD_NAME='경부고속도로'")
변경 전 변경 후

8. Toolbox

8.1 Toolbox 사용하기

  • 사용 가능한 Toolbox 기능 리스트 확인하기
    • QgsApplication.processingRegistry().algorithms() : 사용 가능한 알고리즘 전체 객체 리스트
    • alg.provider().name() : 알고리즘 provider 이름
    • alg.name() : 알고리즘 이름
    • alg.displayName() : 알고리즘 표시 이름
from qgis import processing
for alg in QgsApplication.processingRegistry().algorithms():
  print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))
GDAL:aspect --> 경사 방향
GDAL:assignprojection --> 투영체 적용
GDAL:buffervectors --> 벡터 버퍼
GDAL:buildvirtualraster --> 가상 래스터 생성
GDAL:buildvirtualvector --> 가상 벡터 생성
GDAL:cliprasterbyextent --> 범위로 래스터 자르기
...
  • 알고리즘 이름으로 찾기
[x.id() for x in QgsApplication.processingRegistry().algorithms() if "buffer" in x.id()]
['gdal:buffervectors', 'gdal:onesidebuffer', 'grass7:r.buffer', 'grass7:r.buffer.lowmem', 'grass7:v.buffer', 'native:buffer', 'native:bufferbym', 'native:multiringconstantbuffer', 'native:singlesidedbuffer', 'native:taperedbuffer', 'native:wedgebuffers', 'qgis:variabledistancebuffer', 'saga:fixeddistancebuffer', 'saga:rasterbuffer', 'saga:rasterproximitybuffer', 'saga:thresholdrasterbuffer', 'saga:variabledistancebuffer']
  • 알고리즘 사용방식 확인
processing.algorithmHelp("native:buffer")
버퍼 (native:buffer)

이 알고리즘은 고정 또는 동적 거리를 사용해서 입력 레이어의 모든 객체에 대해 버퍼 영역을 계산합니다.

둥근 오프셋을 생성하는 경우 선분 파라미터가 사분원을 비슷하게그리는 데 사용할 라인 선분의 개수를 제어합니다.

선끝(end cap) 스타일 파라미터는 버퍼 내부에서 라인 끝부분을 어떻게 처리할지 제어합니다.

결합 스타일 파라미터는 오프셋이 라인의 모서리에 적용될 경우 결합 부위를 둥글게(round) 할지, 마이터(miter)로 할지, 비스듬하게(bevel) 할지 지정합니다.

마이터 제한 파라미터는 마이터 결합 스타일에만 적용할 수 있으며, 마이터 결합 부위를 생성하는 경우 사용할 오프셋 곡선으로부터의 최대 거리를 제어합니다.


----------------
Input parameters
----------------

INPUT: 입력 레이어

    Parameter type:    QgsProcessingParameterFeatureSource

    Accepted data types:
        - str: 레이어 ID
        - str: 레이어 이름
        - str: 레이어 원본
        - QgsProcessingFeatureSourceDefinition
        - QgsProperty
        - QgsVectorLayer

DISTANCE: 거리

    Parameter type:    QgsProcessingParameterDistance

    Accepted data types:
        - int
        - float
        - QgsProperty

SEGMENTS: 선분

    반올림한 오프셋을 생성하는 경우, 선분 파라미터가 사분원을 근사치로 계산하기 위해 사용할 라인 선분의 개수를 제어합니다.

    Parameter type:    QgsProcessingParameterNumber

    Accepted data types:
        - int
        - float
        - QgsProperty

END_CAP_STYLE: 선끝 스타일

    Parameter type:    QgsProcessingParameterEnum

    Available values:
        - 0: 둥글게
        - 1: 평평하게
        - 2: 정사각형

    Accepted data types:
        - int
        - str: int를 표현하는 문자열, 예. '1'
        - QgsProperty

JOIN_STYLE: 이음새 스타일

    Parameter type:    QgsProcessingParameterEnum

    Available values:
        - 0: 둥글게
        - 1: 마이터(miter)
        - 2: 비스듬하게(bevel)

    Accepted data types:
        - int
        - str: int를 표현하는 문자열, 예. '1'
        - QgsProperty

MITER_LIMIT: 마이터 제한

    Parameter type:    QgsProcessingParameterNumber

    Accepted data types:
        - int
        - float
        - QgsProperty

DISSOLVE: 결과물 디졸브

    Parameter type:    QgsProcessingParameterBoolean

    Accepted data types:
        - bool
        - int
        - str
        - QgsProperty

OUTPUT: 산출물

    Parameter type:    QgsProcessingParameterFeatureSink

    Accepted data types:
        - str: 대상 벡터 파일, 예. 'd:/test.shp'
        - str: 임시 메모리 레이어에 결과를 저장하는 'memory :'
        - str: 벡터 공급자 ID 접두사 및 대상 URI 사용, 예. PostGIS 테이블에 결과를 저장하는 'postgres:…'
        - QgsProcessingOutputLayerDefinition
        - QgsProperty

----------------
Outputs
----------------

OUTPUT:  <QgsProcessingOutputVectorLayer>
    산출물
  • 알고리즘 사용하기(buffer)
result = processing.run("native:buffer",{'INPUT':'mNode','DISTANCE':10,'SEGMENTS':5, 'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2, 'DISSOLVE':False,'OUTPUT':'memory:'})
QgsProject.instance().addMapLayer(result['OUTPUT'])

8.2 Toolbox 만들기

참고

  • 파이썬을 활용한 지리공간 분석 마스터하기
  • 예제는 국가표준노드링크로 진행
  • QGIS 파이썬 라이브러리는 여기 참고
반응형
반응형

1. 데이터 읽기

  • 벡터 데이터 읽기
    • read_file(*.shp) : shp파일 읽어서 geoDataFrame에 저장
import geopandas as gpd
df = gpd.read_file(r'data/MOCT_LINK.shp',encoding='CP949')
df.head()

LINK_ID F_NODE T_NODE LANES ROAD_RANK ROAD_TYPE ROAD_NO ROAD_NAME ROAD_USE MULTI_LINK CONNECT MAX_SPD REST_VEH REST_W REST_H LENGTH REMARK geometry
0 2630193301 2630076801 2630076901 1 106 000 391 화악산로 0 0 000 60 0 0.0 0 1410.192910 None LINESTRING (245889.208 602540.103, 245884.524 ...
1 2630193001 2630076801 2630076701 1 106 003 391 화악산로 0 0 000 60 0 0.0 0 12.137670 None LINESTRING (245881.843 602537.719, 245885.114 ...
2 2630193101 2630076701 2630076801 1 106 003 391 화악산로 0 0 000 60 0 0.0 0 12.326808 None LINESTRING (245893.460 602528.496, 245889.209 ...
3 2630192801 2630076701 2630076601 1 106 000 391 화악산로 0 0 000 60 0 0.0 0 364.089006 None LINESTRING (245885.688 602526.172, 245886.272 ...
4 2630192901 2630076601 2630076701 1 106 000 391 화악산로 0 0 000 60 0 0.0 0 373.389143 None LINESTRING (246066.042 602242.391, 246069.650 ...
  • Geometry 표시 하기
%matplotlib inline
df.plot(color='black')

2. geoDataFrame 구조 확인

  • type 확인
type(df)
geopandas.geodataframe.GeoDataFrame
  • row, column 크기 확인
df.shape
(1331, 18)
  • column명 확인
df.columns
Index(['LINK_ID', 'F_NODE', 'T_NODE', 'LANES', 'ROAD_RANK', 'ROAD_TYPE',
       'ROAD_NO', 'ROAD_NAME', 'ROAD_USE', 'MULTI_LINK', 'CONNECT', 'MAX_SPD',
       'REST_VEH', 'REST_W', 'REST_H', 'LENGTH', 'REMARK', 'geometry'],
      dtype='object')
  • Geometry type 확인
df.geom_type
0         LineString
1         LineString
2    MultiLineString
3         LineString
4         LineString
dtype: object

3. 좌표계 확인

  • 좌표계 확인
df.crs
<Projected CRS: PROJCS["ITRF2000_Central_Belt_60",GEOGCS["GCS_ITRF ...>
Name: ITRF2000_Central_Belt_60
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Transverse Mercator
Datum: International Terrestrial Reference Frame 2000
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
  • 좌표계 변경
merc = df.to_crs({'init':'epsg:4326'})
merc.plot(color='black')

4. 포맷 변환

  • GeoDataFrame을 json 포맷으로 변환
# 데이터가 커서 sample로 1개 데이터만 json으로 변환
df.head(1).to_json()
'{"type": "FeatureCollection", "features": [{"id": "0", "type": "Feature", "properties": {"CONNECT": "000", "F_NODE": "2630076801", "LANES": 1, "LENGTH": 1410.19291001225, "LINK_ID": "2630193301", "MAX_SPD": 60, "MULTI_LINK": "0", "REMARK": null, "REST_H": 0, "REST_VEH": "0", "REST_W": 0.0, "ROAD_NAME": "\\ud654\\uc545\\uc0b0\\ub85c", "ROAD_NO": "391", "ROAD_RANK": "106", "ROAD_TYPE": "000", "ROAD_USE": "0", "T_NODE": "2630076901"}, "geometry": {"type": "LineString", "coordinates": [[245889.20842293755, 602540.1031623926], [245884.52438196362, 602550.7076639828], [245880.5845997197, 602562.441661735], [245877.37697640058, 602577.5560817202], [245874.80402102915, 602590.923142394], [245870.08166344027, 602608.6555760134], [245866.14860266628, 602619.1390588756], [245861.72407481138, 602627.9941839982], [245850.47753049707, 602649.6933080491], [245843.77653209324, 602663.2882877537], [245840.21593625328, 602674.2739925202], [245836.2741359605, 602686.3831435101], [245833.96002093932, 602698.125878997], [245831.74675705165, 602714.3711385361], [245830.48617608764, 602739.5004281478], [245829.72781690946, 602764.257252305], [245829.78834916753, 602799.523103117], [245829.18865190083, 602818.0280274062], [245827.4648357864, 602836.2767985966], [245825.86203727854, 602855.2765506067], [245823.64944289793, 602871.396758187], [245821.52694263405, 602894.0203128068], [245820.4236764545, 602913.1478050507], [245821.58007512547, 602930.6617291515], [245823.60645547815, 602949.1807706261], [245825.40895076204, 602962.8214614922], [245826.09405003514, 602974.9554847644], [245825.65368618732, 602987.0834576795], [245823.97894215182, 602996.2034733664], [245820.42842643222, 603005.3134050232], [245811.20958754103, 603022.0212222983], [245803.40463366537, 603031.6085023026], [245795.35831640614, 603039.5687689325], [245779.41022924686, 603051.8634824678], [245761.20651050736, 603065.0214546616], [245744.9962197665, 603079.5657484786], [245711.8058337661, 603112.2767922183], [245676.3322525288, 603150.9782006346], [245664.11688441734, 603166.7945190094], [245654.27749311822, 603182.6236104805], [245646.6896994245, 603198.3397514643], [245641.49132830047, 603211.5676374831], [245635.8969615895, 603228.67010031], [245632.192484073, 603243.1565677221], [245627.70139361764, 603264.3917775303], [245625.5082893798, 603276.885489703], [245622.31073984868, 603290.1241318531], [245618.2163164351, 603307.3597129482], [245615.66553742788, 603316.6000712622], [245612.62396502416, 603324.0870205313], [245606.79428443874, 603338.4370064927], [245600.21765677864, 603352.1577013484], [245595.16584919338, 603361.384613071], [245585.18325133767, 603380.5894169775], [245577.72319553196, 603395.8060218558], [245572.39439199885, 603410.0336447023], [245567.3082955352, 603425.6381774854], [245563.35505705074, 603439.873195574], [245557.6094144114, 603461.851987991], [245552.1448012958, 603478.0797602631], [245546.8139796491, 603492.6825364484], [245543.12428970283, 603504.4178705714], [245537.91112339962, 603520.3968842088], [245527.91978245924, 603541.2273534045], [245501.14465575287, 603590.2300162113], [245479.78640583853, 603631.8835551282], [245454.11186099224, 603685.5191653903], [245430.8179966427, 603738.292187912], [245408.38805921804, 603793.195787896], [245399.2431610517, 603819.4081672801]]}}]}'
  • GeoDataFrame을 json 파일로 저장
    • driver : 출력 파일 포맷
df.to_file(driver='GeoJSON', filename='link.geojson')
  • 파일 출력 시 가능한 포맷
    • GeoPandas는 Fiona 라이브러리를 이용해 출력
    • 아래 명령을 통해 Fiona에서 지원하는 모든 driver 확인 가능
import fiona
fiona.supported_drivers

5. 데이터 다루기

  • 첫 번째 데이터 확인
df.loc[0]
LINK_ID                                              2180383101
F_NODE                                               2180132301
T_NODE                                               2180132201
LANES                                                         1
ROAD_RANK                                                   101
ROAD_TYPE                                                   000
ROAD_NO                                                      17
ROAD_NAME                                              서울문산고속도로
ROAD_USE                                                      0
MULTI_LINK                                                    0
CONNECT                                                     101
MAX_SPD                                                      50
REST_VEH                                                      5
REST_W                                                        0
REST_H                                                        0
LENGTH                                                  169.152
REMARK                                                     None
geometry      LINESTRING (310943.9252600061 4166093.43858687...
Name: 0, dtype: object
  • 특정 컬럼 데이터 조회
df['ROAD_NAME']
0       서울문산고속도로
1         개화동로8길
2         개화동로8길
3            당산길
4              -
          ...
1326    서울문산고속도로
1327    서울문산고속도로
1328    서울문산고속도로
1329    서울문산고속도로
1330    서울문산고속도로
Name: ROAD_NAME, Length: 1331, dtype: object
  • 컬럼 필터링
서울문산고속도로 = df[df['ROAD_NAME'] == '서울문산고속도로']
서울문산고속도로
  • Geometry 보기
서울문산고속도로.plot(figsize=(7,7))

6 공간 Join

  • 서울시에 포함되는 Link 찾기 예제
    • 행정경계는 여기에서 다운로드 가능
df2 = gpd.read_file(r'data/TL_SCCO_CTPRVN.shp')
df2.plot()

seoul_link = gpd.sjoin(df, seoul, op='within')

참고

반응형
반응형

1. Shapely로 Geometry 생성하기

  • Polygon 생성하기(튜플)
from shapely.geometry import Polygon
p1 = Polygon(((1,2), (5,3), (5,7), (1,9), (1,2)))
p1

p2 = Polygon(((6,6), (7,6), (10,4), (11,8), (6,6)))
p2

  • Point생성하기(좌표값)
from shapely.geometry import Point
point = Point(2.0, 2.0)
point

  • Line 생성하기
from shapely.geometry import LineString
line = LineString([(0,0), (10,10)])
line

  • Linear Ring 생성하기
from shapely.geometry.polygon import LinearRing
ring = LinearRing([(0,0), (3,3), (3,0)])
ring

  • MultiPoint 생성하기
from shapely.geometry import MultiPoint
points = MultiPoint([(0,0), (3,3)])
points

  • MultiLine 생성하기
from shapely.geometry import MultiLineString
coords = MultiLineString([((0,0), (1,1)),((-1,0), (1,0))])
coords

  • MultiPolygon 생성하기
from shapely.geometry import MultiPolygon
polygons = MultiPolygon([p1,p2])
polygons

2. Shapely 공간 함수


# 면적 구하기
print(p1.area)
# 경계 구하기
print(p1.bounds)
# 길이 구하기
print(p1.length)
# geometry Type 구하기
print(p1.geom_type)
22.0
(1.0, 2.0, 5.0, 9.0)
19.59524158061724
Polygon

3. Shapely로 JSON Geometry 읽기

  • Json 파일 생성 및 shape를 이용해 로딩
import json
from shapely.geometry import mapping, shape
jData = json.loads('{"type":"Polygon","coordinates":[[[1,1],[1,3],[3,3]]]}')
p = shape(jData)
p

  • geometry를 json형태로 변환
mapping(p)
{'type': 'Polygon', 'coordinates': (((1.0, 1.0), (1.0, 3.0), (3.0, 3.0), (1.0, 1.0)),)}

4 Fiona 데이터 읽기

  • shp 파일 정보 표출
print(len(c))
print(c.driver)
print(c.crs)
100
ESRI Shapefile
{'proj': 'tmerc', 'lat_0': 38, 'lon_0': 127, 'k': 1, 'x_0': 200000, 'y_0': 600000, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}
  • shp 파일 데이터 읽기
import fiona
c = fiona.open(r'data/node.shp')
rec = next(iter(c))
print('keys = ', rec.keys())
print('type = ', rec['type'])
print('prop = ', rec['properties'])
print('id   = ', rec['id'])
print('geom = ', rec['geometry'])
keys =  dict_keys(['type', 'id', 'properties', 'geometry'])
type =  Feature
prop =  OrderedDict([('NODE_ID', '1100025100'), ('NODE_TYPE', '101'), ('NODE_NAME', '하계5,6단지앞교차로'), ('TURN_P', '1'), ('REMARK', None)])
id   =  0
geom =  {'type': 'Point', 'coordinates': (205838.0880999937, 559449.3750999967)}

5 데이터 다루기

  • shape file 읽고 첫 번째 Feature 정보 가져오기
import fiona
with fiona.open(r'data/node.shp') as src:
    print(src[0])
{'type': 'Feature', 'id': '0', 'properties': OrderedDict([('NODE_ID', '1100025100'), ('NODE_TYPE', '101'), ('NODE_NAME', '하계5,6단지앞교차로'), ('TURN_P', '1'), ('REMARK', None)]), 'geometry': {'type': 'Point', 'coordinates': (205838.0880999937, 559449.3750999967)}}

참고

반응형
반응형

1. OGR 라이브러리

  • OGR의 3가지 Component
    • OGR 배치 명령 : 벡터 데이터 정보 요약 및 변환
    • 파이썬 스크립트 : 미리 개발해 놓은 GIS 분석 기능
    • OGR 라이브러리 : 파이썬에서 사용할 수 있는 OGR 라이브러리

2. OGR 배치 명령

  • OGR은 벡터데이터 요약 정보 제공, 변환 등을 위한 배치 명령을 제공
  • 대표적으로 아래와 같은 배치 명령 존재
    • ogrinfo : 벡터 데이터 요약정보 확인
    • ogr2ogr : 벡터 파일을 다른 포맷으로 변경
    • ogrtindex : 벡터 타일 생성
    • ogr2vrt : 벡터 타일 생성(orgtindex 보다 광범위하게 사용됨)
  • 터미널에서 명령을 수행함
    • Jupyter notebook에서 명령어 앞에 !를 붙이면 터미널 명령 수행 가능

2.1 ogrinfo

  • 파일 입출력 가능 포맷 조회
!ogrinfo --formats
Supported Formats:
...
  ESRI Shapefile -vector- (rw+v): ESRI Shapefile
  KML -vector- (rw+v): Keyhole Markup Language (KML)
  GeoJSON -vector- (rw+v): GeoJSON
...
  • 파일의 요약 정보 보기
    • data 폴더의 link 파일 요약 정보 보기
    • -so : 요약 정보만 보겠다는 옵션
    • -so가 없으면 shp의 모든 정보를 표시함
!ogrinfo -so "data" MOCT_LINK
INFO: Open of `data'
      using driver `ESRI Shapefile' successful.

Layer name: MOCT_LINK
Metadata:
  DBF_DATE_LAST_UPDATE=2021-03-30
Geometry: Line String
Feature Count: 528232
Extent: (101766.537600, 67516.493200) - (546275.776100, 665746.517600)
Layer SRS WKT:
PROJCS["ITRF2000_Central_Belt_60",
    GEOGCS["GCS_ITRF_2000",
        DATUM["International_Terrestrial_Reference_Frame_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6656"]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",38],
    PARAMETER["central_meridian",127],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",200000],
    PARAMETER["false_northing",600000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
Data axis to CRS axis mapping: 1,2
LINK_ID: String (10.0)
F_NODE: String (10.0)
T_NODE: String (10.0)
LANES: Integer (4.0)
ROAD_RANK: String (3.0)
ROAD_TYPE: String (3.0)
ROAD_NO: String (5.0)
ROAD_NAME: String (30.0)
ROAD_USE: String (1.0)
MULTI_LINK: String (1.0)
CONNECT: String (3.0)
MAX_SPD: Integer (4.0)
REST_VEH: String (3.0)
REST_W: Integer (4.0)
REST_H: Integer (4.0)
LENGTH: Real (18.12)
REMARK: String (30.0)
!ogrinfo "KML_Sample.kml" -summary
INFO: Open of `KML_Sample.kml'
      using driver `KML' successful.
1: Placemarks (3D Point)
2: Highlighted Icon (3D Point)
3: Paths (3D Line String)
4: Google Campus (3D Polygon)
5: Extruded Polygon (3D Polygon)
6: Absolute and Relative (3D Polygon)

2.2 ogr2ogr

  • shp 파일을 GeoJson파일로 변환
!ogr2ogr -f "GeoJson" "output.json" "data/MOCT_LINK.shp"

3. 파이썬 스크립트

  • GDAL 설치 경로에 GIS 분석에 사용할 수 있는 스크립트 제공
    • 경로 예시(anaconda일 때) : C:\Users\82103\Anaconda3\envs\geospatial\Lib\site-packages\GDAL-3.0.2-py3.6-win-amd64.egg-info\scripts
  • scripts 예시
    • ogrmerge.py : shp파일을 merge 하여 다른 파일로 변환함
  • 사용 방법
    • 터미널 또는 Jupyter notebook에서 사용 방식에 따라 명령 수행
#경로 이동
cd "C:\Users\82103\Anaconda3\envs\geospatial\Lib\site-packages\GDAL-3.0.2-py3.6-win-amd64.egg-info\scripts"

#터미널에서 수행 할 때
ogrmerge.py -f GPKG -o merged.gpkg "MOCT_LINK.shp"

# Jupyter Note Book에서 수행할 때
%run ogrmerge.py -f GPKG -o merged.gpkg "MOCT_LINK.shp"
  • 참고 : %run
    • Jupyter notebook의 매직명령
    • 외부 스크립트를 Jupyter notebook에서 실행 할 때 사용

4. OGR 라이브러리

  • OGR 라이브러리는 2개의 주요 모듈로 구성
    • ogr : 벡터 지오메트리를 주로 다룸
    • osr : 투영에 관해 다룸
  • OGR 이 제공하는 7개의 클래스
    • Geometry : 기하학 좌표 연산/변환 제공, Spatial Reference System을 포함
    • Spatial Reference : 투영과 좌표계 관련 기능 정의
    • Feature : Geometry와 속성정보를 가진 Feature를 정의(Polygon 등의 객체)
    • Feature Class Definition : Feature의 스키마 정보를 정의
    • Layer : Feature들을 모아놓은 Layer 정의
    • Dataset : OGRLayer 객체를 하나 이상 포함하는 File 또는 Database을 표현
    • Drivers : 포맷 변환을 정의하며, 모든 Drivers는 GDALDriverManager에서 관리함

4.1 OGR로 객체 생성하기

  • OGR로 Point, Line, Polygon과 같은 벡터 Geometry 생성 가능
  • Polygon 생성 예제
from osgeo import ogr
r = ogr.Geometry(ogr.wkbLinearRing)
r.AddPoint(1,1)
r.AddPoint(5,1)
r.AddPoint(5,5)
r.AddPoint(1,5)
r.AddPoint(1,1)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(r)
print(poly.ExportToWkt())
POLYGON ((1 1 0,5 1 0,5 5 0,1 5 0,1 1 0))

4.2 GeoJSON으로 부터 객체 생성하기

  • GeoJSON 구조로 부터 객체 생성 가능
  • GeoJSON으로 부터 Polygon 생성 예제
from osgeo import ogr
geojson = """{"type":"Polygon","coordinates":[[[1,1],[5,1],[5,5],[1,5],[1,1]]]}"""
polygon = ogr.CreateGeometryFromJson(geojson)
print(polygon)
POLYGON ((1 1,5 1,5 5,1 5,1 1))

4.3 Geometry 연산 수행

  • Polygon Geometry 연산 수행
# 1. 면적 구하기
print("폴리곤 면적은 ", polygon.Area(), "입니다.")

# 2. 중심점 구하기
cen = polygon.Centroid()
print(cen)

# 3. 경계 구하기
boundary = polygon.GetBoundary()
print(boundary)

# 4. Point 가 Polygon 내부에 존재하는지 확인
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10, 10)
polygon.Contains(point)

# 5.  Buffer 구하기
buffer = polygon.Buffer(2)
폴리곤 면적은  16.0 입니다.
POINT (3 3)
LINESTRING (1 1,5 1,5 5,1 5,1 1)
False

4.4 Shape 파일에 Polygon Data 쓰기

  • 아래 순서로 진행
    • 좌표계 설정
    • Shape file 생성
    • Layer 생성
    • feature 생성 및 Geometry 입력
    • Layer에 Feature 입력
import osgeo.ogr, osgeo.osr

# 1. 좌표계 설정
spatialReference = osgeo.osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

# 2. ShapeFile 생성
driver = osgeo.ogr.GetDriverByName('ESRI ShapeFile')
shapeData = driver.CreateDataSource('polygon.shp')

# 3. layer 생성
layer = shapeData.CreateLayer('polygon_layer', spatialReference, osgeo.ogr.wkbPolygon)
layerDefinition = layer.GetLayerDefn()

# 4. Feature 생성 및 Geometry 입력
featureIndex = 0
feature = osgeo.ogr.Feature(layerDefinition)
feature.SetGeometry(polygon)
feature.SetFID(featureIndex)

# 5. Layer에 Feature 입력
layer.CreateFeature(feature)

# 6. 생성 확인
!ogrinfo polygon.shp
INFO: Open of `polygon.shp'
      using driver `ESRI Shapefile' successful.
1: polygon (Polygon)

4.5 feature 다루기

  • Shape 파일을 읽어 특정 필드명 출력
from osgeo import ogr
import os
shapefile = r'data/node.shp'
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile,0)
layer = dataSource.GetLayer()

# 데이터가 많아서 Sample로 5개 필드만 출력
idx = 0
for feature in layer:
    print(idx, ", ", feature.GetField("NODE_NAME"))
    idx += 1
    if(idx == 5):
        break
0 ,  하계5,6단지앞교차로
1 ,  공릉터널남측
2 ,  미성아파트5동앞
3 ,  유토피아빌딩
4 ,  청구빌라4단지401동
  • Shape 파일의 모든 필드 이름 나열하기
from osgeo import ogr
source = ogr.Open(r'data/MOCT_LINK.shp')
layer = source.GetLayer()
schema = []
ldefn = layer.GetLayerDefn()
for n in range(ldefn.GetFieldCount()):
    fdefn = ldefn.GetFieldDefn(n)
    schema.append(fdefn.name)
print(schema)
['LINK_ID', 'F_NODE', 'T_NODE', 'LANES', 'ROAD_RANK', 'ROAD_TYPE', 'ROAD_NO', 'ROAD_NAME', 'ROAD_USE', 'MULTI_LINK', 'CONNECT', 'MAX_SPD', 'REST_VEH', 'REST_W', 'REST_H', 'LENGTH', 'REMARK']
  • Feature 개수 구하기
from osgeo import ogr
import os
shapefile = r'data/MOCT_LINK.shp'
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
featureCount = layer.GetFeatureCount()
print("count = ", featureCount)
count =  528232

4.6 좌표계 다루기

  • CRS 정보 조회하기
    • Layer로부터 구하기
    • Geometry로부터 구하기
# 1) Layer로부터 구하기
spatialRef = layer.GetSpatialRef()
print(spatialRef)
PROJCS["ITRF2000_Central_Belt_60",
    GEOGCS["GCS_ITRF_2000",
        DATUM["International_Terrestrial_Reference_Frame_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6656"]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",38],
    PARAMETER["central_meridian",127],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",200000],
    PARAMETER["false_northing",600000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
# 2) Geometry로부터 구하기
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef2 = geom.GetSpatialReference()
print(spatialRef2)
PROJCS["ITRF2000_Central_Belt_60",
    GEOGCS["GCS_ITRF_2000",
        DATUM["International_Terrestrial_Reference_Frame_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6656"]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",38],
    PARAMETER["central_meridian",127],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",200000],
    PARAMETER["false_northing",600000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

참고

반응형
반응형

1. 벡터 데이터 타입

벡터 데이터 타입 특징
ShapeFile - 가장 많이 쓰는 포맷으로 3개의 파일로 구성
   - shp : Geometry 정보 보유
   - shx : Geomtery를 빠르게 찾을 수 있게 인덱스 보유
   - dbf : 속성 정보 보유
GeoJson - Json 기반의 공간 정보 포맷으로 key, value 쌍으로 구성
- 확장자는 .json 또는 .geojson을 사용
KML - XML을 기반으로 태그 형태의 포맷
- 확장자는 압축된 .kmz 형태로 배포되는 경우가 많음
- WGS84 에 정의된 경도, 위도 좌표계를 사용
GeoPackage - 벡터와 래스터 데이터를 모두 지원하는 개방형 데이터 포맷
- 모바일 사용자를 위한 포맷
- 데이터와 메타데이터 테이블을 결합한 확장 SQLite3 데이터베이스 파일(.gpkg) 사용

2.래스터 데이터 타입

래스터 데이터 타입 특징
ECW - 항공 및 위성 사진용으로 많이 사용
- 압축된 이미지 포맷으로 화질을 유지하면서 압축율이 높음
Esri grid - 래스터 파일에 속성 데이터를 추가한 파일 포맷
GeoTIFF - GIS 및 원격 탐사 App을 위한 표준 이미지 파일
반응형
반응형

1. postgresql 설치(postGIS 포함)

2. 파이썬에서 PostgreSQL 연결하기

  • psycopg2 설치하기
conda install -c anaconda psycopg2

3. 데이터 베이스 연결하기

  • geospatial DB에 연결 후 points 테이블 생성하는 예제
import psycopg2, pprint

connection = psycopg2.connect(database="geospatial", user="postgres", password="******")
cursor = connection.cursor()

cursor.execute("CREATE TABLE points (id SERIAL PRIMARY KEY, name VARCHAR(255), location GEOMETRY)")
connection.commit()
  • psycopg2.connect
    • 데이터베이스명, 사용자ID, PW를 입력하여 DB에 연결
  • connection.cursor
    • 데이터베이스와 통신할 수 있는 cursor 생성
  • cursor.execute
    • 인자로 입력된 SQL문을 실행
  • connection.commit
    • 변경 내용을 DB에 저장

4. 테이블에 데이터 추가

  • WKT형태의 좌표를 Geometry로 저장
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(100 200)'))")             
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(200 400)'))")
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(700 300)'))")
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(400 300)'))")
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(100 400)'))")
cursor.execute("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('POINT(100 300)'))")
connection.commit() 

cursor.execute("SELECT * from points") 
cursor.fetchall() 
  • WKT 대신 shapely를 이용해 geometry를 컨트롤 하는 방법
    • shapely의 Point 객체이용
    • Point 객체는 wkt 함수를 이용해 WKT로 변환 가능
from shapely.geometry import Point, MultiPoint
p = [Point(100,200), Point(100,100)]
for pp in p:
    s = ("INSERT INTO points (name, location) VALUES ('p1',ST_GeomFromText('{}'))").format(pp.wkt)
    cursor.execute(s)
connection.commit() 

cursor.execute("SELECT * from points") 
cursor.fetchall() 
  • 주피터 노트북에서 이미지 띄우기
    • 좌표 리스트들이 표시됨
MultiPoint(p)

5. 테이블 조회

  • points를 조회 하는 예제
    • fetchall() 을 이용하면 select 된 모든 데이터 가져옴
    • fetchone을 이용하면 1개의 데이터만 가져옴
cursor.execute("SELECT * from points") 
data=cursor.fetchall() 
data
[(1, 'p1', '010100000000000000000059400000000000006940'),
 (2, 'p1', '010100000000000000000069400000000000007940'),
 (3, 'p1', '01010000000000000000E085400000000000C07240'),
 (4, 'p1', '010100000000000000000079400000000000C07240'),
 (5, 'p1', '010100000000000000000059400000000000007940'),
 (6, 'p1', '010100000000000000000059400000000000C07240'),
 (8, 'p1', '010100000000000000000059400000000000006940'),
 (9, 'p1', '010100000000000000000059400000000000005940')]
  • geomtery가 WKB 형태로 출력
    • shapely를 이용해 WKB-> Point -> WKT 변환 하기
from shapely.wkb import loads
aPoint=loads(data[0][2],hex=True) 
aPoint.wkt
'POINT (100 200)'

6. 좌표계 설정/변경하기

cursor.execute("SELECT UpdateGeometrySRID('points','location',4326)") 
cursor.execute("SELECT Find_SRID('public','points','location')") 
cursor.fetchall()

cursor.execute("SELECT code, ST_AsTexT(ST_Transform(location,3857)) from points") 
cursor.fetchone()
  • updateGeomterySRID
    • Geometry에 좌표계 설정함
  • Find_SRID
    • Geomtery가 어떤 좌표계 인지 보여줌
  • ST_Transform
    • 좌표계를 변경함

7. KNN을 이용한 지정한 점에 가장 가까운 3개 점 얻기

cursor.execute("SELECT name, ST_Distance(location::geography, ST_GeometryFromText('POINT(-106.591838300225 35.1555000000615)')::geography) as d from points ORDER BY location<->ST_GeometryFromText('POINT(-106.591838300225 35.1555000000615)') LIMIT 3") 
cursor.fetchall() 
[('p1', 7120405.75751067),
 ('p1', 16921336.66007583),
 ('p1', 16921336.66007583)]

8. line 데이터 다루기

  • line 생성 예제
from shapely.geometry import LineString 
from shapely.geometry import MultiLineString
connection = psycopg2.connect(database="geospatial",user="postgres",password="****")
cursor = connection.cursor() 
cursor.execute("CREATE TABLE lines (id SERIAL PRIMARY KEY, location GEOMETRY)") 
thelines=[] 
thelines.append(LineString([(-106.635585,35.086972),(-106.621294,35.124997)])) 
thelines.append(LineString([(-106.498309,35.140108),(-106.497010,35.069488)])) 
thelines.append(LineString([(-106.663878,35.106459),(-106.586506,35.103979)]))
mls=MultiLineString([((-106.635585,35.086972),(-106.621294,35.124997)),((-106.498309,35.140108),(-106.497010,35.069488)),((-106.663878 ,35.106459),(-106.586506,35.103979))])
for a in thelines:    
    cursor.execute("INSERT INTO lines (location) VALUES (ST_GeomFromText('{}'))".format(a.wkt)) 
    connection.commit() 
cursor.execute("SELECT id, ST_AsTexT(location) from lines") 
data=cursor.fetchall() 
data
connection.commit()
  • 두 선 교차 여부 확인 하기
    • ST_Intersects 함수 사용
cursor.execute("SELECT ST_Intersects(l.location::geography,ll.location::geometry) FROM lines l, lines ll WHERE l.id=1 AND ll.id=3") 
cursor.fetchall()
  • 두 선을 교차하는 점 확인 하기
cursor.execute("SELECT ST_AsText(ST_Intersection(l.location::geography, ll.location::geometry)) FROM lines l, lines ll WHERE l.id=1 AND ll.id=3") 
cursor.fetchall()

9.Polygon 데이터 다루기

  • polygon 생성 예제

    from shapely.geometry import Polygon
    connection = psycopg2.connect(database="pythonspatial",user="postgres", password="postgres") 
    cursor = conectionn.cursor() 
    cursor.execute("CREATE TABLE poly (id SERIAL PRIMARY KEY, location GEOMETRY)") 
    a=Polygon([(-106.936763,35.958191),(-106.944385,35.239293),(-106.452396,35.281908),(-106.407844,35.948708)]) 
    cursor.execute("INSERT INTO poly (location) VALUES (ST_GeomFromText('{}'))".format(a.wkt)) 
    connection.commit()
  • 포인트가 폴리곤 내부에 있는지 확인

    • ST_Intersects 또는 ST_Contains 사용

# 참고

반응형
반응형

1. GDAL/OGR

  • GDAL 과 OGR 이란?
    • GDAL(Geospatial Data Abstraction Library : 래스터 데이터 처리에 많이 사용하는 라이브러리
    • OGR : 벡터 포맷 데이터 처리에 많이 사용하는 라이브러리
  • GDAL 설치 하기
    • 아나콘다 내비게이터를 이용해 설치
      • 가상 환경 생성
      • 생성된 가상환경에 Not installed 설정 -> gdal 검색하여 체크 -> Apply
    • 아나콘다 프롬프트를 이용해 설치
      • 가상 환경 생성
conda create -n geospatial
activate geospatial
conda install gdal

2. GEOS(Geometry Engine Open Source)

  • GEOS란?
    • JTS(Java Topology Suite)의 서브셋과 선택된 기능의 C/C++ 포트
    • JTS의 완전한 기능을 C++로 포함하는 것을 목표
    • QGIS, PostGIS 등 여러 어플리케이션에서 사용됨
    • GDAL과 함께 컴파일해서 OGR의 모든 기능을 제공
  • JTS란?
    • Java로 작성된 지리공간 지오메트리 라이브러리
  • GEOS 설치하기
conda install -c conda-forge geos

3. Shapely

  • Shapely란?
    • 평면 형상의 처리와 분석을 위한 파이썬 패키지
    • GEOS 라이브러리 및 JTS 포트의 기능을 사용함
    • 지오메트리 분석만 다룸(공간 파일을 읽고 쓸 수 있는 기능 없음)
    • 8가지 기본 지오메트리 타입 지원
      • 포인트, 멀티 포인트, 멀티 라인 스트링, 멀티 폴리곤, 폴리곤, 다원형, 지오메트리 컬렉션
    • OGR보다 더 파이썬에 가깝고 직관적인 인터페이스, 최적화가 잘 되어 있음
  • Shapely 설치하기
conda install -c scitools shapely

4. Fiona

  • Fiona란?
    • OGR API로 데이터 포맷을 읽고 쓰는데 사용
    • OGR 보다 파이썬에 가깝고 신뢰성이 높고 오류 발생률이 낮아 OGR대신 사용
    • 벡터 데이터에 관한 공간 정보를 표현하기 위해 WKT와 WKB라는 두 개의 마크업 언어 사용
    • Shapely와 같은 다른 파이썬 라이브러리와 잘 결합됨
    • 벡터 데이터 복사 시 메모리를 많이 사용(C의 포인터가 아니라 파이썬의 오브젝트 사용하기 때문)
  • Fiona 설치
conda install -c conda-forge fiona

5. pyshp

  • pyshp(Python Shapefile Library) 란?
    • Shape File을 다루기 위한 파이썬 라이브러리
    • Shape File만 다룬다면 GDAL을 쓰는 것 보다 pyshp을 사용하는 것을 추천
  • pyshp 설치
conda install pyshp

6. pyproj

  • pyproj 란?
    • 도표 변환과 측지 연산을 수행하는 파이썬 패키지
    • PROJ.4 함수에 파이썬 인터페이스를 제공하기 위한 Cython 래퍼로 파이썬에서 기존 C코드 라이브러리에 직접 접근 가능
  • PROJ.4 란?
    • 여러 좌표계 중에서 데이터를 변환하는 투영 라이브러리
    • GDAL과 OGR에서도 이용 가능
    • 인기 있는 이유
      • 다양한 좌표계를 지원
      • Rasterio, GeoPandas 라이브러리 모두 pyproj를 사용하고, 그에 따라 PROJ.4 기능을 사용함
  • pyproj 설치
conda install -c conda-forge pyproj

7. Rasterio

  • Rasterio 란?
    • 래스터 데이터를 위한 GDAL 및 Numpy 기반 파이썬 라이브러리
    • Web과 App을 위한 온라인 지도 제공업체 맵박스 위성팀에서 제공하는 오픈소스
  • Rasterio 설치
conda config --add channels conda-forge
conda install rasterio

8. GeoPandas

  • GeoPandas 란?
    • 벡터 데이터 작업을 위한 파이썬 라이브러리
    • pandas와 유사하고 공간 데이터 처리에 유용하게 사용됨
  • GeoPandas 설치
conda install -c conda-forge geopandas

# 참고

  • 파이썬을 활용한 지리 공간 분석 마스터하기
반응형

+ Recent posts

반응형