Proyeksi Ortho menghasilkan artefak

Saya mencoba membuat tampilan seperti bola menggunakan qgis dan proyeksi "dunia dari luar angkasa" http://spatialreference.org/ref/sr-org/6980/ (terutama proyeksi-ortho). ArcGIS membungkus bentuk dengan benar tetapi QGIS (2.01) menghasilkan artefak jahat.

Saya harus menghasilkan bola secara teratur dengan sudut yang berbeda, jadi apakah ada orang di luar sana yang punya ide bagaimana memperbaiki masalah ini?

Seperti yang dikatakan Andre, agar ini berfungsi, Anda harus memotong layer Anda sebelum memproyeksikannya. Andre menjelaskan metode manual , yang bekerja dengan baik untuk banyak kasus: Proyeksikan shapefile Anda ke proyeksi sama rata dengan parameter yang sama dengan proyeksi ortografis, buat lingkaran kliping yang menutupi belahan bumi yang akan terlihat dalam proyeksi ortografi, dan jepitkan shapefile dengan itu. Namun, metode itu memerlukan sedikit upaya manual, dan tidak berfungsi untuk semua parameter proyeksi, karena memproyeksikan ke proyeksi azimuth yang sama dapat menyebabkan masalah yang sama seperti memproyeksikan ke proyeksi ortografis.

Berikut ini skrip (sekarang juga tersedia sebagai plugin QGIS Klip ke Belahan ) yang menggunakan pendekatan yang sedikit berbeda: Lapisan kliping dibuat dalam sistem referensi koordinat dari shapefile asli dengan memproyeksikan lingkaran dari ortografis ke CRS sumber, tetapi juga tambahan pastikan untuk menutupi seluruh belahan yang terlihat, termasuk tiang yang terlihat.

Inilah yang terlihat seperti lapisan kliping untuk proyeksi ortografis yang berpusat pada 30 ° N, 110 ° E:

Script kemudian klip lapisan yang dipilih saat ini dengan lapisan kliping, dan menambahkan lapisan yang dihasilkan ke proyek. Lapisan itu kemudian dapat diproyeksikan ke proyeksi ortografis, baik dengan cepat atau dengan menyimpannya dalam CRS ortografis:

Ini skripnya. Pastikan untuk menyimpannya di jalur Python Anda, misalnya sebagai 'cliportho.py'. Kemudian Anda dapat mengimpornya menggunakan konsol QGIS Python import cliportho. Untuk memotong sebuah layer, panggil cliportho.doClip(iface, lat=30, lon=110, filename='A.shp').

import numpy as np
from qgis.core import *
import qgis.utils

import sys, os, imp


def doClip(iface, lat=30, lon=110, filename='result.shp'):
    sourceLayer = iface.activeLayer()

    sourceCrs = sourceLayer.dataProvider().crs()

    targetProjString = "+proj=ortho +lat_0=" + str(lat) + " +lon_0=" + str(lon) + "+x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
    targetCrs = QgsCoordinateReferenceSystem()
    targetCrs.createFromProj4(targetProjString)

    transformTargetToSrc = QgsCoordinateTransform(targetCrs, sourceCrs).transform

    def circlePolygon(nPoints=20, radius=6370000, center=[0,0]):
        clipdisc = QgsVectorLayer("Polygon?crs=epsg:4326", "Clip disc", "memory")
        angles = np.linspace(0, 2*np.pi, nPoints, endpoint=False)
        circlePoints = np.array([ transformTargetToSrc(QgsPoint(center[0]+np.cos(angle)*radius, center[1]+np.sin(angle)*radius)) for angle in angles ])
        sortIdx = np.argsort(circlePoints[:,0])
        circlePoints = circlePoints[sortIdx,:]
        circlePoints = [ QgsPoint(point[0], point[1]) for point in circlePoints ]
        circlePoints.extend([QgsPoint(180,circlePoints[-1][1]), QgsPoint(180,np.sign(lat)*90), QgsPoint(-180,np.sign(lat)*90), QgsPoint(-180,circlePoints[0][1])])
        circle = QgsFeature()
        circle.setGeometry(QgsGeometry.fromPolygon( [circlePoints] ) )
        clipdisc.dataProvider().addFeatures([circle])
        QgsMapLayerRegistry.instance().addMapLayer(clipdisc)
        return clipdisc

    auxDisc = circlePolygon(nPoints = 3600)

    ###### The clipping stuff
    ## Code taken from the fTools plugin

    vproviderA = sourceLayer.dataProvider()
    vproviderB = auxDisc.dataProvider()

    inFeatA = QgsFeature()
    inFeatB = QgsFeature()
    outFeat = QgsFeature()

    fitA = vproviderA.getFeatures()

    nElement = 0  
    writer = QgsVectorFileWriter( filename, 'UTF8', vproviderA.fields(),
                                  vproviderA.geometryType(), vproviderA.crs() )

    index = QgsSpatialIndex()
    feat = QgsFeature()
    index = QgsSpatialIndex()
    fit = vproviderB.getFeatures()
    while fit.nextFeature( feat ):
        index.insertFeature( feat )

    while fitA.nextFeature( inFeatA ):
      nElement += 1
      geom = QgsGeometry( inFeatA.geometry() )
      atMap = inFeatA.attributes()
      intersects = index.intersects( geom.boundingBox() )
      first = True
      found = False
      if len( intersects ) > 0:
        for id in intersects:
          vproviderB.getFeatures( QgsFeatureRequest().setFilterFid( int( id ) ) ).nextFeature( inFeatB )
          tmpGeom = QgsGeometry( inFeatB.geometry() )
          if tmpGeom.intersects( geom ):
            found = True
            if first:
              outFeat.setGeometry( QgsGeometry( tmpGeom ) )
              first = False
            else:
              try:
                cur_geom = QgsGeometry( outFeat.geometry() )
                new_geom = QgsGeometry( cur_geom.combine( tmpGeom ) )
                outFeat.setGeometry( QgsGeometry( new_geom ) )
              except:
                GEOS_EXCEPT = False
                break
        if found:
          try:
            cur_geom = QgsGeometry( outFeat.geometry() )
            new_geom = QgsGeometry( geom.intersection( cur_geom ) )
            if new_geom.wkbType() == 0:
              int_com = QgsGeometry( geom.combine( cur_geom ) )
              int_sym = QgsGeometry( geom.symDifference( cur_geom ) )
              new_geom = QgsGeometry( int_com.difference( int_sym ) )
            try:
              outFeat.setGeometry( new_geom )
              outFeat.setAttributes( atMap )
              writer.addFeature( outFeat )
            except:
              FEAT_EXCEPT = False
              continue
          except:
            GEOS_EXCEPT = False
            continue
    del writer

    resultLayer = QgsVectorLayer(filename, sourceLayer.name() + " - Ortho: Lat " + str(lat) + ", Lon " + str(lon), "ogr")
    QgsMapLayerRegistry.instance().addMapLayer(resultLayer)

Anda harus memotong data poligon Anda ke belahan dunia yang terlihat, karena QGIS tidak melakukannya sendiri.

Saya menulis tutorial di sini:

Ke mana perginya poligon setelah memproyeksikan peta di QGIS?

EDIT

Gambar yang Anda perlihatkan sebenarnya bukan proyeksi ortho, karena menunjukkan seluruh dunia, dan tidak hanya separuh yang terlihat seperti yang terlihat dari luar angkasa. Untuk peta dunia, pemotongannya sedikit lebih mudah, seperti dijelaskan di sini:

QGIS menampilkan file bentuk negara dunia yang berpusat di laut pasifik menggunakan Robinson, Miller Cylindrical atau proyeksi lainnya