1. Rasterio로 래스터 데이터 읽고 쓰기

• Rasterio로 파일 읽기

import rasterio
dataset = rasterio.open(r'Natural_Earth_quick_start\50m_raster\NE1_50M_SR_W\NE1_50M_SR_W.tif')

• 밴드 수 확인

dataset.count

3

• 열 수(width) 조회

dataset.width

10800

• 높이(height) 조회

dataset.height

5400

• 공간 경계(boundbox) 조회

dataset.bounds

BoundingBox(left=-179.99999999999997, bottom=-89.99999999998201, right=179.99999999996405, top=90.0)

• CRS 조회

dataset.crs

CRS.from_epsg(4326)

• 밴드1 데이터 조회

band1 = dataset.read(1)

• 이미지 보기

%matplotlib inline
from matplotlib import pyplot
pyplot.imshow(dataset.read(1))
pyplot.show()

2. GDAL로 래스터 데이터 읽고 쓰기

2.1 터미널에서 GDAL 명령 수행하기

• 터미널 명령어 이며, 상기 내용 처럼 !를 포함하면 JupyterNoteBook에서 명령 가능하다.
• 지원되는 모든 파일 형식 조회

!gdalinfo --formats

Supported Formats:
  VRT -raster- (rw+v): Virtual Raster
...

• 파일 정보 요약(CRS 포함)

!gdalinfo "Natural_Earth_quick_start\50m_raster\NE1_50M_SR_W\NE1_50M_SR_W.tif"

Driver: GTiff/GeoTIFF
Files: Natural_Earth_quick_start\50m_raster\NE1_50M_SR_W\NE1_50M_SR_W.tif
Size is 10800, 5400
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-179.999999999999972,90.000000000000000)
Pixel Size = (0.033333333333330,-0.033333333333330)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2014:10:18 09:32:38
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CC 2014 (Macintosh)
  TIFFTAG_XRESOLUTION=342.85699
  TIFFTAG_YRESOLUTION=342.85699
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Band 1 Block=10800x1 Type=Byte, ColorInterp=Red
Band 2 Block=10800x1 Type=Byte, ColorInterp=Green
Band 3 Block=10800x1 Type=Byte, ColorInterp=Blue

• GeoTiff를 JPEG로 변환

!gdal_translate -of JPEG "Natural_Earth_quick_start\50m_raster\NE1_50M_SR_W\NE1_50M_SR_W.tif" NE1_50M_SR_W.jpg

Input file size is 10800, 5400
0...10...20...30...40...50...60...70...80...90...100 - done.

• GDAL을 사용해 GeoPackage 파일 정보 보기
  • Sample 파일 다운로드

!gdalinfo gdal_sample_v1.2_no_extensions.gpkg

Driver: GPKG/GeoPackage
Files: gdal_sample_v1.2_no_extensions.gpkg
Size is 512, 512
Subdatasets:
  SUBDATASET_1_NAME=GPKG:gdal_sample_v1.2_no_extensions.gpkg:byte_png
  SUBDATASET_1_DESC=byte_png - byte_png
  SUBDATASET_2_NAME=GPKG:gdal_sample_v1.2_no_extensions.gpkg:byte_jpeg
  SUBDATASET_2_DESC=byte_jpeg - byte_jpeg
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

2.2 GDAL 라이브러리 사용

2.2.1 GDAL 읽기/쓰기/투영 확인/변경

• sample 데이터에서 New Mexico Color Shaded Relief 데이터 다운
• 데이터 읽기

from osgeo import gdal
nmtif = gdal.Open(r'Data/nm_relief_color_tif/nm_relief_color.tif')
print(nmtif.GetMetadata())

{'AREA_OR_POINT': 'Area', 'TIFFTAG_DATETIME': '2002:12:18  8:10:06', 'TIFFTAG_RESOLUTIONUNIT': '2 (pixels/inch)', 'TIFFTAG_SOFTWARE': 'IMAGINE TIFF Support\nCopyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile: etif.c $ $Revision: 1.9.3.3 $ $Date: 2002/07/29 15:51:11EDT $', 'TIFFTAG_XRESOLUTION': '96', 'TIFFTAG_YRESOLUTION': '96'}

• 투영 정보 확인

nmtif.GetProjection()

'PROJCS["NAD83 / UTM zone 13N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26913"]]'

• 투영 변경 하기

from osgeo import osr
p = osr.SpatialReference()
p.ImportFromEPSG(26913)
nmtif.SetProjection(p.ExportToWkt())

• 새 TIF 파일로 저장

geoTiffDriver = "GTiff"
driver=gdal.GetDriverByName(geoTiffDriver)
out = driver.CreateCopy('copy.tif', nmtif, strict = 0)

2.2.2 데이터 다루기

• 밴드 개수 구하기

nmtif.RasterCount

3

• 1번 밴드 데이터 보기

band = nmtif.GetRasterBand(1)
values = band.ReadAsArray()
values

array([[254, 255, 255, ..., 250, 249, 247],
       [255, 255, 255, ..., 251, 250, 249],
       [255, 255, 255, ..., 252, 250, 249],
       ...,
       [255, 255, 255, ..., 255, 255, 254],
       [255, 255, 255, ..., 254, 253, 253],
       [255, 255, 255, ..., 252, 251, 251]], dtype=uint8)

values[1100,1100]

216

• 특정 픽셀 컬러값 구하기

one = nmtif.GetRasterBand(1).ReadAsArray()
two = nmtif.GetRasterBand(2).ReadAsArray()
three = nmtif.GetRasterBand(3).ReadAsArray()
print(str(one[1100,1100]) + "," + str(two[1100,1100]) + "," + str(three[1100,1100]))

216,189,157

• 특정 band 평균, 표준편차 얻기

one = nmtif.GetRasterBand(1)
one.ComputeBandStats()

(225.05771967375847, 34.08382839593031)

• 특정 band 최대/최소값 구하기

print("최대값 = ", str(one.GetMaximum()))
print("최소값 = ", str(one.GetMinimum()))

최대값 =  255.0
최소값 =  0.0

• 특정 band의 description 확인/변경

print("변경 전 = ", str(one.GetDescription()))
one.SetDescription("The Red Band")
print("변경 후 = ", str(one.GetDescription()))

변경 전 =  Band_1
변경 후 =  The Red Band

• Jupyter Notebook에서 래스터 파일 보기

import numpy as np
from matplotlib.pyplot import imshow
%matplotlib inline
data_array=nmtif.ReadAsArray()
x=np.array(data_array[0])
image = x.reshape(x.shape[0], x.shape[1])
imshow(image, cmap='gist_earth')

2.2.3 GDAL을 이용한 래스터 생성

• Sample Data 생성
  • numpy array 이용

raster = np.array([[10,10,1,10,10,10,10],
                  [1,1,1,50,10,10,50],
                 [10,1,1,51,10,10,50],
                 [1,1,1,1,50,10,50]])
raster

array([[10, 10,  1, 10, 10, 10, 10],
       [ 1,  1,  1, 50, 10, 10, 50],
       [10,  1,  1, 51, 10, 10, 50],
       [ 1,  1,  1,  1, 50, 10, 50]])

• 래스터 파일 속성 설정

# 좌하단 모서리 좌표 설정
coord=(-106,629773, 35.105389)

# 폭 , 높이 설정
w, h = 10, 10

# 파일명 설정
name = "new.tif"

• 래스터 파일에 쓰기

# 1. driver 가져오기
d=gdal.GetDriverByName("GTiff")

# 2. file 설정
rows = raster.shape[0]
cols = raster.shape[1]
numOfBands = 1
output = d.Create(name, cols, rows, numOfBands, gdal.GDT_UInt16)
output.SetGeoTransform((coord[0], w, 0, coord[1], 0 ,h))

# 3. 데이터 쓰기
output.GetRasterBand(1).WriteArray(raster)

# 4. 좌표계 설정
outsr = osr.SpatialReference()
outsr.ImportFromEPSG(4326)
output.SetProjection(outsr.ExportToWkt())

# 5. 파일에 쓰기
output.FlushCache()

• 결과 확인

data = output.ReadAsArray()
w,h = 4,7
image = data.reshape(w,h)
imshow(image, cmap='Blues')
data

array([[10, 10,  1, 10, 10, 10, 10],
       [ 1,  1,  1, 50, 10, 10, 50],
       [10,  1,  1, 51, 10, 10, 50],
       [ 1,  1,  1,  1, 50, 10, 50]], dtype=uint16)

# 참고

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

1. 벡터 데이터 타입

2.래스터 데이터 타입

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

# 참고

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