Orthoprojektion erzeugt Artefakte

I`m versucht , einen kugelähnliche Ansicht mit qgis zu schaffen und die „Welt aus dem Weltraum“ -Projektion http://spatialreference.org/ref/sr-org/6980/ (essentialy eine ortho-Projektion). ArcGIS bricht die Formen korrekt um, QGIS (2.01) erzeugt jedoch unangenehme Artefakte.

Ich muss die Globen regelmäßig mit verschiedenen Winkeln herstellen. Hat jemand eine Idee, wie man dieses Problem beheben kann?

Wie Andre sagte, damit dies funktioniert, müssen Sie Ihre Ebene beschneiden, bevor Sie sie projizieren. Andre beschreibt eine manuelle Methode , die in vielen Fällen gut funktioniert: Projizieren Sie Ihr Shapefile auf eine azimutale äquidistante Projektion mit denselben Parametern wie die orthografische Projektion, erstellen Sie einen Beschneidungskreis, der die Halbkugel bedeckt, die in der orthografischen Projektion sichtbar sein wird, und Clip das Shapefile damit. Diese Methode erfordert jedoch ein wenig manuellen Aufwand und funktioniert nicht bei allen Projektionsparametern, da die Projektion auf eine azimutale Projektion mit gleichem Abstand zu ähnlichen Problemen führen kann wie die Projektion auf eine orthografische Projektion.

Hier ist ein Skript (jetzt auch als Clip-to-Hemisphere-QGIS-Plug-in verfügbar ), das einen etwas anderen Ansatz verfolgt: Im Koordinatenreferenzsystem des ursprünglichen Shapefiles wird eine Clipping-Ebene erstellt, indem ein Kreis vom orthografischen zum Quell-CRS projiziert wird Stellen Sie sicher, dass die gesamte sichtbare Hemisphäre, einschließlich des sichtbaren Pfostens, bedeckt ist.

So sieht die Clipping-Ebene für eine orthografische Projektion aus, die auf 30 ° N und 110 ° E zentriert ist:

Das Skript schneidet dann die aktuell ausgewählte Ebene mit der Schnittebene ab und fügt die resultierende Ebene dem Projekt hinzu. Diese Ebene kann dann direkt oder durch Speichern in der orthografischen CRS auf die orthografische Projektion projiziert werden:

Hier ist das Drehbuch. Stellen Sie sicher, dass Sie es in Ihrem Python-Pfad speichern, zum Beispiel als 'cliportho.py'. Anschließend können Sie es mit in die QGIS Python-Konsole importieren import cliportho. Um eine Ebene auszuschneiden, rufen Sie auf 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)

Sie müssen Ihre Polygondaten auf die sichtbare Hälfte des Globus zuschneiden, da QGIS dies nicht alleine tut.

Ich habe hier ein Tutorial geschrieben:

Wohin gingen die Polygone, nachdem sie eine Karte in QGIS projiziert hatten?

BEARBEITEN

Das Bild, das Sie zeigen, ist eigentlich keine Orthoprojektion, da es die ganze Welt zeigt und nicht nur die sichtbare Hälfte aus dem Weltraum. Für Weltkarten ist das Schneiden etwas einfacher, wie hier beschrieben:

QGIS zeigt mit Robinson, Miller Cylindrical oder einer anderen Projektion weltweite Formdateien an, die auf den Pazifischen Ozean ausgerichtet sind