소개

QGIS에 내장된 Python 콘솔을 사용하여 GIS 데이터를 OBJ 파일로 내보냅니다.

※OBJ 파일에 대해서는 Wikipedia 링크 를 제시해 둡니다. 텍스트로 내보낼 수 있는 매우 간단한 포맷이므로, 독자적으로 export할 경우에 취급하기 쉽습니다.

자르기

먼저 대상 범위를 자릅니다. 전회의 데이터를 그대로 사용하면, 과연 거대해져서, 리얼타임 렌더링에서는 대응할 수 없게 되므로, 적당한 범위에 잘라 봐 봅니다.
이번에는 롯폰기 힐즈 주변을 잘라 보겠습니다.
먼저 메뉴에서 [레이어]-[레이어 만들기]-[새 임시 스크래치 레이어]를 선택하고 대화 상자에서 적절한 이름과 형상 유형으로 다각형을 선택합니다.
추가한 새 레이어가 편집 모드에 있는지 확인하고 메뉴에서 편집 > 직사각형 추가 > 영역 범위 직사각형 추가를 수행합니다. 정점과 대각선의 정점을 왼쪽 클릭으로 선택하고 오른쪽 클릭으로 추가 종료가 됩니다.



그림은 알기 쉽게 구형의 묘화를 반투명하게 하고 있습니다.

그런 다음 '벡터' - '공간 연산 도구' - '자르기'를 선택하고 입력에 건물 데이터의 레이어, 오버레이에 지금 만든 사각형 영역이있는 레이어를 선택하고 실행합니다.

이제 내보낼 영역이 잘립니다.



그런 다음이 직사각형 영역의 무게 중심을 찾습니다.
「벡터」-「지오메트리 툴」-「중심」입니다.
가능한 포인트를 피처 정보 표시 도구로 클릭하면 무게 중심의 위도 경도를 알 수 있습니다.

QGIS 매크로? 파이썬 콘솔

이제 Python 콘솔을 사용하여 이 데이터를 3D OBJ 파일로 내보냅니다.
Python 콘솔은 메뉴에서 '플러그인'-'Python 콘솔'로 호출할 수 있습니다.

좌표 변환

현재 좌표계는 WGS84로 (아무것도 바꾸지 않으면 아마) 작업하고 있습니다. 좌표계에 대해 세세한 곳을 설명하기 시작하면 엄청난 양이므로 여기에서는 할애합니다.
WGS84에서 좌표는 위도 경도로 표시됩니다. 여기에서 XYZ의 3차원 직교 좌표계로 변환합니다. 다음의 메소드는, 세계측지계와 좌표변환 라고 하는 책의 프로그램을 참고로 해 작성한 좌표 변환의 메소드입니다. WGS84에서 지구의 표현 방법에 따른 최대 정확한 직교 좌표 변환이 되어 있습니다.
덧붙여 이 계산에서는, 적어도 double의 정밀도로 실시할 필요가 있습니다. Python의 부동 소수점은 double이어야하므로 별로 신경 쓰지 않아도 될까 생각합니다만, Unity의 Vector3 클래스 등은 float로 수치를 가지고 있으므로, 이 프로그램을 이식할 때는 주의가 필요합니다.
종종 지리 좌표를 다룰 때는 double로 다루는 것이 더 이상 실패하지 않습니다.

def BLH2XYZ(b,l,h):
  # WGS84
  a = 6378137.0
  f = 1.0 / 298.257223563

  b = math.radians(b)
  l = math.radians(l)

  e2 = f * (2.0 - f)
  N = a / math.sqrt(1.0 - e2 * math.pow(math.sin(b), 2.0))

  X = (N + h) * math.cos(b) * math.cos(l);
  Y = (N + h) * math.cos(b) * math.sin(l);
  Z = (N * (1 - e2) + h) * math.sin(b);

  return (X,Y,Z)

OBJ 내보내기

그리고 위의 좌표 변환 메소드를 포함한 형태로 OBJ 내보내기 프로그램을 작성합니다.

def BLH2XYZ(b,l,h):
  # WGS84
  a = 6378137.0
  f = 1.0 / 298.257223563

  b = math.radians(b)
  l = math.radians(l)

  e2 = f * (2.0 - f)
  N = a / math.sqrt(1.0 - e2 * math.pow(math.sin(b), 2.0))

  X = (N + h) * math.cos(b) * math.cos(l);
  Y = (N + h) * math.cos(b) * math.sin(l);
  Z = (N * (1 - e2) + h) * math.sin(b);

  return (X,Y,Z)

layer = iface.activeLayer()
features = layer.getFeatures()

i = 0
j = 0
fp = open("test.obj",mode='w')
fp.write("g "+str(i)+"\n")

# ※前に取得しておいた矩形領域の重心に変える
cx = 139.72957
cy = 35.66021

ox,oy,oz = BLH2XYZ(cy,cx,0)

for feature in features:
    # print(feature.id())
    mp = feature.geometry().asMultiPolygon()
    try:
        height = int(feature['height'])*2
        if height < 1:
            height = 2
    except:
        height = 5

    for polygon in mp:
        for polyline in polygon:
            i=i+1
            prev_p_index = j

            for point in polyline:
                x,y,z = BLH2XYZ(point.y(),point.x(),0)
                x2,y2,z2 = BLH2XYZ(point.y(),point.x(),height)
                x = x - ox
                y = y - oy
                z = z - oz
                x2 = x2 - ox
                y2 = y2 - oy
                z2 = z2 - oz

                s = math.radians(-cx)
                rx = x * math.cos(s) - y * math.sin(s)
                ry = x * math.sin(s) + y * math.cos(s)

                s = math.radians(-cy)
                rxx = rx * math.cos(s) - z * math.sin(s)
                rz  = rx * math.sin(s) + z * math.cos(s)

                s = math.radians(-cx)
                rx2 = x2 * math.cos(s) - y2 * math.sin(s)
                ry2 = x2 * math.sin(s) + y2 * math.cos(s)

                s = math.radians(-cy)
                rxx2 = rx2 * math.cos(s) - z2 * math.sin(s)
                rz2  = rx2 * math.sin(s) + z2 * math.cos(s)


                fp.write("v "+str(ry)+" "+str(rxx)+" "+str(rz)+"\n")
                fp.write("v "+str(ry2)+" "+str(rxx2)+" "+str(rz2)+"\n")
                j=j+1

            current = j - prev_p_index
            for num in range(current):
                p1 = (2*num+1)
                p2 = (2*num+2)
                p3 = (2*num+3)
                p4 = (2*num+4)

                if p3 > current*2:
                    p3 = p3 - current*2

                if p4 > current*2:
                    p4 = p4 - current*2

                p1 = p1 + prev_p_index*2
                p2 = p2 + prev_p_index*2
                p3 = p3 + prev_p_index*2
                p4 = p4 + prev_p_index*2

                fp.write("f "+str(p1)+" "+str(p2)+" "+str(p3)+" "+str(p1)+"\n")
                fp.write("f "+str(p4)+" "+str(p3)+" "+str(p2)+" "+str(p4)+"\n")

fp.close()
print("Done")

이 프로그램 내에서, 처리로서는 방금전의 지리 좌표로부터 직교 좌표로의 좌표 변환과, 원점의 이동과 전체의 회전, OBJ에의 기입을 실시하고 있습니다.

결과