Ortho Projection produce artefactos

Estoy tratando de crear una vista en forma de esfera usando qgis y la proyección "mundo desde el espacio" http://spatialreference.org/ref/sr-org/6980/ (esencialmente una orto-proyección). ArcGIS envuelve las formas correctamente, pero QGIS (2.01) produce artefactos desagradables.

Tengo que producir los globos de forma regular con diferentes ángulos, ¿alguien tiene alguna idea de cómo solucionar este problema?

Como dijo Andre, para que esto funcione, deberá recortar su capa antes de proyectarla. Andre describe un método manual , que funciona bien para muchos casos: proyecte su archivo de forma en una proyección equidistante azimutal con los mismos parámetros que la proyección ortográfica, cree un círculo de recorte que cubra el hemisferio que será visible en la proyección ortográfica, y recorta el archivo de forma con eso. Sin embargo, ese método requiere un poco de esfuerzo manual y no funciona para todos los parámetros de proyección, ya que proyectar en una proyección equidistante acimutal puede conducir a problemas similares a los de proyectar en una proyección ortográfica.

Aquí hay un script (ahora también disponible como el complemento Clip to Hemisphere QGIS ) que adopta un enfoque ligeramente diferente: se crea una capa de recorte en el sistema de referencia de coordenadas del archivo de forma original al proyectar un círculo desde el CRS ortográfico al origen, pero adicionalmente asegurándose de cubrir todo el hemisferio visible, incluido el polo visible.

Así es como se ve la capa de recorte para una proyección ortográfica centrada en 30 ° N, 110 ° E:

Luego, el script recorta la capa seleccionada actualmente con la capa de recorte y agrega la capa resultante al proyecto. Esa capa se puede proyectar en la proyección ortográfica, ya sea sobre la marcha o guardándola en el CRS ortográfico:

Aquí está el guión. Asegúrese de guardarlo en su ruta de Python, por ejemplo como 'cliportho.py'. Luego puede importarlo en la consola QGIS Python usando import cliportho. Para recortar una capa, llame 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)

Debe recortar los datos de su polígono a la mitad visible del globo, porque QGIS no lo hace por sí solo.

Escribí un tutorial aquí:

¿A dónde fueron los polígonos después de proyectar un mapa en QGIS?

EDITAR

La imagen que muestra en realidad no es una proyección ortográfica, ya que muestra todo el mundo, y no solo la mitad visible vista desde el espacio exterior. Para los mapas mundiales, el corte es un poco más fácil, como se describe aquí:

QGIS muestra archivos de formas de países del mundo centrados en el océano Pacífico usando Robinson, Miller Cylindrical u otra proyección