Орто-Проекция производит артефакты

Я пытаюсь создать сферическое представление, используя qgis и проекцию «мир из космоса» http://spatialreference.org/ref/sr-org/6980/ (по существу, ортопроекция). ArcGIS правильно оборачивает фигуры, но QGIS (2.01) создает неприятные артефакты.

Я должен регулярно производить глобусы под разными углами, поэтому у кого-нибудь есть идеи, как решить эту проблему?

Как сказал Андре, чтобы это сработало, вам нужно обрезать слой перед его проецированием. Андре описывает ручной метод , который хорошо работает во многих случаях: спроецируйте ваш шейп-файл на азимутальную эквидистантную проекцию с теми же параметрами, что и в ортографической проекции, создайте обтравочный круг, охватывающий полусферу, которая будет видна в ортографической проекции, и обрежьте шейп-файл с этим. Однако этот метод требует значительных ручных усилий и не работает для всех параметров проекции, поскольку проецирование на азимутальную эквидистантную проекцию может привести к тем же проблемам, что и проекция на ортографическую проекцию.

Вот скрипт (теперь также доступный как плагин Clip to Hemisphere QGIS ), который использует немного другой подход: слой отсечения создается в системе координат координат исходного шейп-файла, проецируя окружность от ортогонального к исходному CRS, но дополнительно следя за тем, чтобы покрыть все видимое полушарие, включая видимый полюс.

Вот так выглядит отсечный слой для ортографической проекции с центром в 30 ° с.ш., 110 ° в.д .:

Затем сценарий обрезает текущий выбранный слой с помощью слоя обрезки и добавляет полученный слой в проект. Затем этот слой может быть спроецирован на ортографическую проекцию, либо на лету, либо сохранением его в ортографическом CRS:

Вот сценарий. Обязательно сохраните его в своем пути Python, например, как «cliportho.py». Затем вы можете импортировать его в консоли QGIS Python, используя import cliportho. Чтобы обрезать слой, позвоните 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)

Вы должны обрезать данные полигона до видимой половины земного шара, потому что QGIS не делает это само по себе.

Я написал учебник здесь:

Куда делись полигоны после проецирования карты в QGIS?

РЕДАКТИРОВАТЬ

Картинка, которую вы показываете, на самом деле не является ортопедической, так как показывает весь мир, а не только видимую половину, видимую из космоса. Для карт мира разрезание немного проще, как описано здесь:

QGIS отображает файлы формы мировой страны, центрированные на Тихом океане, используя проекцию Робинсона, Миллера, Цилиндрическую или другую проекцию