Projekcja orto produkuje artefakty

Usiłuję stworzyć widok podobny do kuli za pomocą qgis i projektu „world from space” http://spatialreference.org/ref/sr-org/6980/ (w zasadzie projekcja orto). ArcGIS poprawnie otacza kształty, ale QGIS (2.01) wytwarza nieprzyjemne artefakty.

Muszę regularnie produkować globusy pod różnymi kątami, więc czy ktoś ma pomysł, jak rozwiązać ten problem?

Jak powiedział Andre, aby to zadziałało, musisz przyciąć warstwę przed jej wyświetleniem. Andre opisuje metodę ręczną , która sprawdza się w wielu przypadkach: rzutuj swój plik kształtu na azymutalną równoodległą projekcję o takich samych parametrach jak projekcja ortograficzna, utwórz okrąg obcinający, który obejmie półkulę, która będzie widoczna w projekcji ortograficznej, i przypnij do niej plik kształtu. Jednak ta metoda wymaga sporo wysiłku ręcznego i nie działa dla wszystkich parametrów projekcji, ponieważ projekcja na azymutalną równoodległą projekcję może prowadzić do podobnych problemów jak projekcja do projekcji ortograficznej.

Oto skrypt (teraz dostępny również jako wtyczka Clip to Hemisphere QGIS ), który przyjmuje nieco inne podejście: Warstwa przycinająca jest tworzona w układzie odniesienia do współrzędnych oryginalnego pliku kształtu poprzez rzutowanie okręgu z ortografii do źródłowego CRS, ale dodatkowo upewniając się, że zakryła całą widoczną półkulę, w tym widoczny biegun.

Tak wygląda warstwa obcinania dla rzutu ortograficznego wyśrodkowanego na 30 ° N, 110 ° E:

Skrypt następnie przycina aktualnie wybraną warstwę warstwą przycinającą i dodaje warstwę wynikową do projektu. Tę warstwę można następnie rzutować na rzut ortograficzny w locie lub zapisując w ortograficznym CRS:

Oto skrypt. Pamiętaj, aby zapisać go w ścieżce Pythona, na przykład jako „cliportho.py”. Następnie możesz zaimportować go do konsoli QGIS Python import cliportho. Aby przyciąć warstwę, zadzwoń 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)

Musisz przyciąć swoje dane wielokąta do widocznej połowy globu, ponieważ QGIS nie robi tego sama.

Napisałem tutorial tutaj:

Gdzie poszły wielokąty po projekcji mapy w QGIS?

EDYTOWAĆ

Pokazany obraz nie jest w rzeczywistości projekcją orto, ponieważ pokazuje cały świat, a nie tylko widzialną połowę widzianą z kosmosu. W przypadku map świata cięcie jest nieco łatwiejsze, jak opisano tutaj:

QGIS wyświetla pliki kształtów świata na środku na Oceanie Spokojnym za pomocą Robinsona, Millera Cylindrycznego lub innej projekcji