GISPython

Mapping Shapefile on Raster Map

, , , ,
  • raster map 위에 polygon을 생성하는 과정을 정리했습니다.
  • 간단한 작업이지만 4 GB 용량의 raster file을 불러오는 과정과 좌표를 맞추는 과정이 중요합니다.

1. geotiff 읽기

  • geotiff는 지형과 건물 등 여러 정보를 담은 raster image입니다.
  • 건물 그림자가 담긴 geotiff 파일을 읽고 그 위에 polygon을 그리고자 합니다.

1.1. matplotlib

  • geotiff는 이미지파일이기 때문에 imagej에서 읽을 수도 있습니다.

  • matplotlib에서도 이미지를 다룰 수 있습니다.

  • plt.imread()기능을 이용해 읽을 수 있지만 문제가 있습니다.

  • 너무 큰 그림은 읽어들이지 못합니다.

  • 심지어 DOS 공격이라고까지 의심합니다.

1.2. gdal

wikipedia:GeoTIFF
GDAL
GDAL install on ubuntu

  • GDAL은 래스터와 벡터 데이터를 상호 교환하는 translator library입니다.

    • ArcGIS, PostGIS를 비롯해 상당히 많은 프로그램에서 사용하는 de-facto 표준입니다.
    • 그러나 왜인지 설치가 만만치 않습니다.

  • 전에는 geopandas설치 과정에서 잘 됐던 것 같은데 새 서버에서 말썽입니다.

    • 윈도에서는 이런 순서로 설치하면 된다고 하고,
    • 리눅스에서는 이렇게 하라고 합니다.
    • 하루 오전을 희생한 끝에 겨우 설치했습니다.

  • 파일을 엽니다.

1
2
3
4
from osgeo import gdal
gdal.UseExceptions()

ds = gdal.Open("파일명")

2. raster data 접근

GDAL: Raster API Tutorial

  • Raster API로 파일의 데이터에 접근할 수 있습니다.
  • .GetDriver()로 파일의 형태를 알 수 있습니다.
1
print(f"Driver: {ds.GetDriver().ShortName}/{ds.GetDriver().LongName}")

Driver: GTiff/GeoTIFF

  • .RasterXSize와 .RasterYSize로 그림의 크기를 알 수 있습니다.
  • .RasterCount로 파일에 담긴 raster band의 수를 알 수 있습니다.
1
print(f"Size is {ds.RasterXSize} x {ds.RasterYSize} x {ds.RasterCount}")

Size is 27001 x 36000 x 1

  • .GetProjection()으로 사용된 좌표계를 확인할 수 있습니다.
1
print(f"Projection is {ds.GetProjection()}")

Projection is PROJCS[“Korea 2000 / Central Belt 2010”,GEOGCS[“Korea 2000”,DATUM[“Geocentric_datum_of_Korea”,SPHEROID[“GRS 1980”,6378137,298.257222101,AUTHORITY[“EPSG”,“7019”]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[“EPSG”,“6737”]],PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,“8901”]],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,“9122”]],AUTHORITY[“EPSG”,“4737”]],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”]],AUTHORITY[“EPSG”,“5186”]]

3. raster data 가져오기

  • .GetRasterBand()를 사용해 raster data를 가져옵니다.
  • 제 데이터에서는 .RasterCount에서 band가 1개라고 했으니 괄호 안에 1을 넣어줘야 합니다.
  • 분명 파이썬 코드이지만 GDAL은 숫자를 1부터 셉니다.
1
band = ds.GetRasterBand(1)
  • band에 저장된 raster band data를 numpy array로 가져옵니다.
  • 제가 익숙한 방식이고, 다른 데이터와 결합해야 하기 때문입니다.
  • 가져온 데이터의 차원을 출력해 잘 읽혔는지 확인합니다.
1
2
3
irrad = band.ReadAsArray()
nrows, ncols = irrad.shape
print(nrows, ncols)

36000 27001

4. 좌표 매핑

matplotlib: origin and extent in imshow

  • numpy array의 index가 0이라고 Korea Belt 2000의 좌표도 0은 아닙니다.

  • 2차원 배열의 좌표와 실제 지점의 좌표를 매핑해야 합니다.

  • .GetGeoTransform() 명령으로 변환 인자들을 추출합니다.

  • raster band에서 가져오는 것이 아닙니다.

1
2
3
4
5
x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
print(x0, dx, dxdy, y0, dydx, dy)

x1 = x0 + dx * ncols  # x limit
y1 = y0 + dy * nrows  # y limit

222100.0 1.0 0.0 434400.0 0.0 -1.0

  • 이제 raster data를 화면에 그려봅시다.
  • matplotlib의 extent는 가로세로축에 실제 데이터값을 찍어줍니다.
1
2
3
4
fig, ax = plt.subplots(figsize=(20, 20))
ax.imshow(irrad, cmap="gist_gray", extent=[x0, x1, y1, y0])
ax.set_aspect("equal", "box")
plt.show()

  • 위 정보를 이용해 position을 index로 변환하는 함수를 만들 수 있습니다.
1
2
3
4
5
6
7
# position to index
def pos2idx(X, Y, ds=ds):
    x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
    X_idx = int((X-x0)/dx)  # x index of X coordinate
    Y_idx = int((Y-y0)/dy)  # y index of Y coordinate
    
    return X_idx, Y_idx

5. Polygon 얹기

Create Shapefile from Polygon
geopandas.read_file

  • 살펴볼 지점을 확대해서 그립니다.
  • 건물들 데이터가 담긴 데이터로부터 목표 건물의 좌표를 추출하고,
  • 이 건물을 중심으로 가로세로 100 pixel 범위의 그림을 그립니다.
1
2
3
4
5
6
7
idx_range = 100 #pixel

fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(irrad[Y_idx-idx_range: Y_idx+idx_range, X_idx-idx_range:X_idx+idx_range], 
          cmap="gist_gray", 
          extent=[idx_x-idx_range, idx_x+idx_range, idx_y-idx_range, idx_y+idx_range])
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from shapely.geometry import mapping, Polygon
import fiona

shppath = "파일 저장 경로"

def create_mor(x, y, width, height, angle, filename, shppath=shppath):
    h_half = height/2
    w_half = width/2
    a_rad = np.deg2rad(angle)
    
    def get_coord(h_half, w_half, a_rad):
        X = x + h_half * np.sin(a_rad) - w_half * np.cos(a_rad)
        Y = y + h_half * np.cos(a_rad) + w_half * np.sin(a_rad)
        return (X, Y)
    
    points = [get_coord(h_half, w_half, a_rad), 
              get_coord(-h_half, w_half, a_rad), 
              get_coord(-h_half, -w_half, a_rad), 
              get_coord(h_half, -w_half, a_rad), 
              get_coord(h_half, w_half, a_rad), 
              ]
    
    schema = {"geometry": "Polygon", "properties": {"id": "int"}}
    
    savefile = os.path.join(shppath, filename)
    with fiona.open(f"{savefile}", "w", "ESRI Shapefile", schema) as c:
        c.write({
            "geometry": mapping(Polygon(points)),
            "properties": {"id": 1}
        })
  • 목표 건물 데이터를 읽고 사각형을 형성합니다.
  • geopandas로 그린 건물은 실제 지도에 들어갈 좌표를 갖고 있습니다.
1
2
3
4
5
6
7
8
9
10
import geopandas as gpd

idx_width = df_dir["OMBB_WIDTH"].loc[idx]
idx_height = df_dir["OMBB_HEIGHT"].loc[idx]
idx_angle = df_dir["OMBB_ANGLE"].loc[idx]

create_mor(idx_x, idx_y, idx_width, idx_height, idx_angle, "idx_100")

idx100 = gpd.read_file('./data/200917_shadow/idx_100/idx_100.shp')
idx100.plot()

  • 이렇게 도형을 shapefile로 만들어주고 geopandas로 읽어옵니다.
  • 예쁘게 겹쳐져 있습니다.


도움이 되셨나요? 카페인을 투입하시면 다음 포스팅으로 변환됩니다

Share