Rasterio 파이썬 패키지를 사용하여 공간 래스터 데이터 처리

이 문서는 지형 공간 래스터 데이터와 관련된 일반적인 작업에 Python 패키지 Rasterio를 사용하는 방법에 대한 빠른 소개를 제공하기 위한 것입니다 . 이것은 주로 나 자신을 파악하는 데 너무 오래 걸린 것들의 모음이므로 관심이 있는 사람과 공유하고 싶습니다.

1. 기본 사항 몇 가지

모든 종류의 지도는 어떤 식으로든 (일반적으로 평면인) 표면에 투영되는 기능 (예: 도로, 건물, 다양한 토지 유형, 지표수 등을 나타냄)으로 구성됩니다. 이러한 기능은 일부 좌표 참조 시스템 ( crs ) 에서 제공되는 좌표를 기반으로 벡터 기능 (또는 일반적으로 모양 , 예를 들어 도로를 설명하는 선, 건물 윤곽을 설명하는 다각형) 으로 나타낼 수 있습니다 .

왼쪽, 아래쪽, 오른쪽 및 위쪽에 대한 제한을 정의하는 일련의 범위 내에서 (잘못 그려진) 래스터 맵을 고려하십시오. 바운드 좌표는 좌표 참조 시스템(crs)으로 정의됩니다. crs에 대한 래스터 맵의 방향과 크기는 "회전 행렬"과 래스터 맵의 원점 좌표의 두 부분으로 구성된 아핀 변환에 의해 정의됩니다(자세한 내용은 아래 참조).

래스터 데이터는 픽셀 배열 형태의 지도 그림으로 생각할 수 있습니다. 모든 픽셀은 아핀(affin) 변환 및 기본 crs를기반으로 지리 공간적 위치에 매핑됩니다각 픽셀의 값은 예를 들어 특정 파장 범위에서 픽셀이 커버하는 영역의 밝기(일부 원격 감지 기기에서 제공)와 같은 일부 척도에 해당합니다. 주어진 crs에서 래스터 맵의 가장자리는래스터 맵의 경계 로 정의됩니다.

래스터 데이터의 표현, 경계, 기본 변환 및 crs는 밀접하게 상호 연관되어 있습니다. 일반적으로 이러한 매개변수 중 하나 이상을 변경하려는 경우 다른 매개변수에도 영향을 미칩니다.

Rasterio는 이 프로세스를 도와줍니다.

2. 데이터 세트

GeoTiff 는 래스터 데이터의 일반적인 이미지 형식입니다. 편리한 GeoTiff 파일이 없는 경우 여기에서 다운로드할 수 있습니다 . 다음 예제는 이 파일을 기반으로 합니다. 다음 코드를 사용하여 GeoTiff 파일을 rasterio의 기본 데이터 구조인 데이터 세트 로 읽을 수 있습니다 .

import rasterio
dataset = rasterio.open('sample.tif')

data다음과 같은 몇 가지 유용한 정보를 전달하는 데이터 세트 개체라고 합니다.

• 기본 좌표 참조 시스템 :

dataset.crs 
CRS.from_epsg(32631)

EPSG:32631은 유럽 일부 지역에서 사용되는 로컬 참조 시스템입니다. http://epsg.io는 사용 가능한 대부분의 좌표 참조 시스템에 대한 정보를 제공합니다. Crs 좌표는 일반적으로 x 및 y로 레이블이 지정됩니다. EPSG:4326( WGS84 와 동일 )의 경우 좌표는 일반적인 경도 및 위도 정의에서 도 단위로 제공됩니다.

• 이 데이터 세트에 포함된 래스터 데이터의 범위 :

dataset.bounds
BoundingBox(left=590520.0, bottom=5780620.0, right=600530.0, top=5790630.0)

제공된 숫자는 단위를 데이터 세트의 crs로 사용합니다. 이 특정 사례( EPSG:32631 )에서 숫자는 crs로 정의된 일부 원점을 기준으로 미터 단위로 지도의 경계 상자를 정의합니다. 다른 crs 정의는 예를 들어 도와 같은 다른 단위를 사용할 수 있습니다.

• 현재 래스터 맵 표현에 사용된 아핀 변환 :

dataset.transform 
Affine(10.0, 0.0, 590520.0, 
       0.0, -10.0, 5790630.0)

Affine 객체 로 구현된 이 변환은 방향(행 또는 열)에서 1픽셀의 변경이 다음과 같은 6개의 매개변수를 사용하여 crs 좌표 변경으로 변환하는 방법을 정의합니다(이 순서대로 모두 변환의 crs와 동일한 단위). 픽셀 열의 변화에 ​​따른 x의 변화(+1픽셀에 대해 +10m), 픽셀 행의 변화에 ​​따른 x의 변화(이 경우 0), x 좌표 원점(590520.0m) , 픽셀 열(이 경우 0)의 변화에 ​​따른 y의 변화, 픽셀 행의 변화에 ​​따른 y의 변화(+1픽셀의 경우 -10m) 및 y 좌표 원점(5790630.0 m 여기). 변환은 "회전 행렬"과 데이터 세트의 crs 좌표 원점(이 문서 상단의 스케치 참조)의 두 부분으로 구성되는 것으로 간주할 수 있습니다.

• 이 래스터 맵에 사용할 수 있는 밴드 또는 채널 수 :

dataset.indexes
(1, 2, 3)

(이 경우 3개의 다른 밴드를 사용할 수 있습니다. 각 밴드는 특정 파장 영역에 대한 그레이스케일 맵을 나타냅니다. Rasterio는 어떤 이유로 Python과 달리 1부터 밴드를 세기 시작합니다.)

다음을 사용하여 데이터 세트의 특정 밴드를 배열로 읽을 수 있습니다.

img = dataset.read(i)

여기서 는 i에서 제공하는 인덱스 중 하나입니다 data.indexes. img실제로는 단지 numpy 배열입니다. 사용 가능한 모든 밴드 배열은 동일한 crs, 변환 및 경계를 공유하기 때문에 동일한 모양(따라서 방향 및 크기 조정)을 가집니다.

예시 이미지를 플로팅해 보겠습니다. 배열 값을 0에서 1까지의 범위로 정규화하는 rasterio.plot메서드 와 matplotlib 축을 생성하고 표시하는 메서드를 사용할 수 있습니다 . BGR 채널을 보유하고 있기 때문에 역 인덱스(또는 채널) 순서를 제공 하지만 RGB를 플롯하려고 합니다.adjust_bandshowimgdatadataset

import numpy as np
from rasterio.plot import show, adjust_band
imgdata = np.array([adjust_band(dataset.read(i)) for i in (3,2,1)])
show(imgdata*3)  # factor 3 to increase brightness

RGB 색상을 사용하는 sample.tif GeoTiff 파일의 플롯.

3. 이미지 좌표와 crs 좌표 간의 변환

지리 공간 데이터를 다루기 때문에 이미지 좌표를 crs 좌표로 또는 그 반대로 변환할 수 있기를 원합니다.

단일 픽셀의 경우 다음을 사용하여 이미지 좌표에서 crs 좌표로 이동할 수 있습니다.

x, y = dataset.xy(row, column)

여기서 x및 데이터세트의 밴드 중 하나에 위치한 픽셀 중심의 crs 좌표를 참조합니다( 정수 y값 으로 제한되지 않음) . 또 다른 방법은 이미지 좌표를 데이터 세트의 실제 변환과 간단히 곱하는 것입니다.rowcolumnrowcolumn

x, y = dataset.transform * (row, column)

두 방법 모두 약간 다른 결과를 제공합니다. 첫 번째 방법은 픽셀 중심 좌표를 제공하는 반면 두 번째 방법은 변환의 원점 좌표에 상대적인 좌표를 제공합니다.

마지막으로 다음을 사용하여 crs 좌표에서 이미지 좌표로 이동할 수 있습니다.

row, column = dataset.index(x, y)

4. 서로 다른 좌표계 간의 좌표 변환

일부 좌표계( 로 정의됨)에 좌표 세트가 있고 dataset.crs이를 새 좌표계로 변환하려는 경우:

from rasterio.warp import transform
from rasterio.crs import CRS
new_crs = CRS.from_epsg(4326)   # standard WGS84 coordinates
new_coo = transform(dataset.crs, new_crs, 
                    xs=[590530], ys=[5790650])
new_coo
([4.326399074221181], [52.25878168834648])

5. 서로 다른 좌표계 간의 벡터 변환

벡터 기능은 매끈한 개체 로 정의할 수 있습니다 . edgelen일부 좌표를 중심으로 x가장자리 길이가 있는 정사각형 영역을 정의하는 데 관심이 있다고 가정해 보겠습니다 y. 관심 영역( roi)을 매끈한 Polygon개체로 정의합니다.

from shapely.geometry import Polygon
x, y, edgelen = 597667, 5787216, 1000  # based on EPSG:32631
roi = Polygon([(x - int(edgelen / 2), y + int(edgelen / 2)),
               (x + int(edgelen / 2), y + int(edgelen / 2)),
               (x + int(edgelen / 2), y - int(edgelen / 2)),
               (x - int(edgelen / 2), y - int(edgelen / 2))])
print(roi)
POLYGON ((590030 5791150, 591030 5791150, 591030 5790150, 590030 5790150, 590030 5791150))

도형을 인쇄하면 도형의 Well-known-text 표현을 사용하여 표시됩니다 .

이제 다음을 기반으로 정의 roi(즉, x좌표 y및 ) 를 가정 하지만 다각형을 다음으로 변환하고 싶습니다 .edgelencrs1crs2

from rasterio.crs import CRS
from rasterio.warp import transform_geom
crs1 = CRS.from_epsg(32631)
crs2 = CRS.from_epsg(4326)
roi2 = transform_geom(crs1, crs2, roi)
roi2
{'type': 'Polygon',
 'coordinates': [[(4.319208814194372, 52.26335783407589),
   (4.333857524152045, 52.26319323507531),
   (4.3335878545153586, 52.254205111266934),
   (4.31894210402294, 52.254369657266636),
   (4.319208814194372, 52.26335783407589)]]}

6. 래스터 데이터 자르기

래스터 데이터를 자르려면 rasterio.mask.mask제거하려는 이미지 영역을 마스크한 다음 제거합니다(키워드 인수를 설정한 경우 crop=True).

from rasterio.mask import mask
crop_img, crop_transform = mask(dataset, shapes=[roi], crop=True)

shape여기에서 "GeoJSON과 같은 딕셔너리 또는 Python 지리 인터페이스 프로토콜을 구현하는 객체"( 문서 참조 )를 나타내며 관심 있는 영역을 정의합니다. roi여기서 다시 다각형을 간단히 사용할 수 있습니다.

이미지에서 잘린 부분을 살펴보겠습니다( crop_img[::-1]샘플 파일인 BGR에서 원래 채널 순서를 RGB로 바꾸는 데 사용함).

crop_img = adjust_band(crop_img)
show(crop_img[::-1])

잘린 이미지는 다채로운 튤립 필드를 보여줍니다.

7. 래스터 데이터 리샘플링

예를 들어 서로 다른 이미징 채널의 해상도를 일치시키기 위해 래스터 데이터를 다시 샘플링해야 하는 경우가 있습니다. 리샘플링은 파일에서 사용할 수 있는 래스터 데이터에서 가장 잘 수행되므로 메모리에 있는 모든 데이터를 먼저 파일에 쓰고 다음 코드를 적용하는 것이 좋습니다(rasterio 문서에서 채택됨 ) . :

from rasterio.enums import Resampling
resample_factor = 0.5 # >1: upsample, <1: downsample

with rasterio.open('sample.tif') as dataset:

    imgdata = dataset.read(out_shape=(dataset.count,
                   int(dataset.height * resample_factor),
                   int(dataset.width * resample_factor)),
                resampling=Resampling.bilinear)

    transform = dataset.transform * dataset.transform.scale(
        (dataset.width / imgdata.shape[-1]),
        (dataset.height / imgdata.shape[-2]))

이 예제 코드는 데이터 세트의 모든 대역을 리샘플링하고 새 transform. 자세한 내용은 관련 래스테리오 문서를 참조하십시오 . 다음은 리샘플링된 이미지입니다.

이 이미지는 원본 sample.tif 와 매우 유사해 보이지만 축 눈금과 레이블을 보면 해당 범위가 원본 이미지 범위의 절반에 불과하다는 것을 알 수 있습니다.

8. 래스터 데이터 재투영

재투영은 다른 crs(예: 다른 좌표계), 다른 변환(예: 래스터 회전) 또는 다른 범위를 기반으로 래스터 데이터에 대한 새로운 표현을 찾는 것을 의미합니다. 여기에는 회전 및 왜곡이 포함될 수 있으므로 재투영이 다소 복잡해집니다.

다음 코드는 new_crs기본 변환을 변경하여 주어진 데이터 세트를 새로운 crs( )로 간단하게 재투영합니다. 코드는 데이터 세트에서 읽고, 새 변환, 래스터 높이 및 너비를 도출하고, 재투영된 데이터( )에 대한 배열을 만들고 new_imgdata마지막으로 재투영을 수행합니다.

from rasterio.warp import calculate_default_transform, reproject
new_crs = CRS.from_epsg(4326)  # WGS84
with rasterio.open('sample.tif') as dataset:
    new_transform, width, height = calculate_default_transform(
        dataset.crs, new_crs, 
        dataset.width, dataset.height, *dataset.bounds)
    imgdata = np.array([dataset.read(i) for i in dataset.indexes])
    new_imgdata = np.zeros(imgdata.shape)
    
    reproject(source=imgdata,
              destination=new_imgdata,
              src_transform=dataset.transform,
              src_crs=dataset.crs,
              dst_transform=new_transform,
              dst_crs=new_crs,
              resampling=Resampling.nearest)

이것은 재투영된 sample.tif 파일입니다.

이것은 WGS84 좌표 참조 시스템으로 재투영된 sample.tif입니다. 이미지는 세로 축을 따라 압축되고 약간 회전되며, 모두 다른 crs에 재투영된 결과입니다.

 

Raserio 참고문서

https://rasterio.readthedocs.io/en/latest/quickstart.html

Python Quickstart — rasterio 1.4dev documentation

https://rasterio.readthedocs.io/en/latest/index.html

Rasterio: access to geospatial raster data — rasterio 1.4dev documentation