Ortho Projection produit des artefacts

J'essaie de créer une vue semblable à une sphère en utilisant qgis et la projection "world from space" http://spatialreference.org/ref/sr-org/6980/ (essentiellement une ortho-projection). ArcGIS enveloppe correctement les formes mais QGIS (2.01) produit des artefacts désagréables.

Je dois produire les globes sur une base régulière avec des angles différents. Est-ce que quelqu'un a une idée de comment résoudre ce problème?

Comme Andre l'a dit, pour que cela fonctionne, vous devrez recadrer votre calque avant de le projeter. Andre décrit une méthode manuelle , qui fonctionne bien dans de nombreux cas: projetez votre fichier de formes sur une projection équidistante azimutale avec les mêmes paramètres que la projection orthographique, créez un cercle de découpage qui couvre l'hémisphère qui sera visible dans la projection orthographique, et découpez le fichier de formes avec cela. Cependant, cette méthode nécessite un peu d'effort manuel et ne fonctionne pas pour tous les paramètres de projection, car la projection sur une projection équidistante azimutale peut entraîner des problèmes similaires à la projection sur une projection orthographique.

Voici un script (désormais également disponible sous la forme du plug-in Clip to Hemisphere QGIS ) qui adopte une approche légèrement différente: Une couche de découpage est créée dans le système de référence de coordonnées du fichier de formes d'origine en projetant un cercle de l'orthographe au CRS source, mais en plus en veillant à couvrir tout l'hémisphère visible, y compris le pôle visible.

Voici à quoi ressemble la couche d'écrêtage pour une projection orthographique centrée sur 30 ° N, 110 ° E:

Le script découpe ensuite le calque actuellement sélectionné avec le calque d’écrêtage et ajoute le calque résultant au projet. Cette couche peut ensuite être projetée sur la projection orthographique, soit à la volée, soit en l'enregistrant dans le CRS orthographique:

Voici le script. Assurez-vous de l'enregistrer dans votre chemin Python, par exemple sous le nom «cliportho.py». Vous pouvez ensuite l'importer dans la console QGIS Python à l'aide de import cliportho. Pour découper un calque, appelez 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)

Vous devez recadrer vos données de polygone dans la moitié visible du globe, car QGIS ne le fait pas seul.

J'ai écrit un tutoriel ici:

Où sont passés les polygones après avoir projeté une carte dans QGIS?

ÉDITER

L'image que vous montrez n'est en fait pas une projection ortho, car elle montre le monde entier, et pas seulement la moitié visible vue de l'espace. Pour les cartes du monde, la découpe est un peu plus facile, comme décrit ici:

QGIS affiche des fichiers de forme de pays du monde centrés sur l'océan Pacifique en utilisant une projection Robinson, Miller cylindrique ou autre