La proyección orto produce artefactos

13

Estoy tratando de crear una vista similar a una esfera usando qgis y la proyección "world from space" enlace (esencialmente una orto-proyección). ArcGIS ajusta las formas correctamente, pero QGIS (2.01) produce artefactos desagradables.

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

    
pregunta user1523709 22.11.2013 - 16:08

2 respuestas

22

Como dijo Andre, para que esto funcione, deberás recortar tu capa antes de proyectarla. Andre describe un método manual , que funciona bien en muchos casos: proyecte su shapefile a 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 recorte el shapefile con ese. Sin embargo, ese método requiere un poco de esfuerzo manual, y no funciona para todos los parámetros de proyección, ya que la proyección a una proyección equidistante azimutal puede conducir a problemas similares a los de la proyección ortográfica.

Aquí hay un script (ahora también disponible como Clip to Hemisphere QGIS plugin ) que tiene un enfoque ligeramente diferente: Se crea una capa de recorte en el sistema de referencia de coordenadas del shapefile original al proyectar un círculo desde el CRS ortográfico al CRS de origen, pero además asegurarse de cubrir todo el hemisferio visible, incluido el polo visible.

Esto es lo que parece la capa de recorte para una proyección ortográfica centrada en 30 ° N, 110 ° E:

Lasecuenciadecomandosluegorecortalacapaseleccionadaactualmenteconlacapaderecorteyagregalacaparesultantealproyecto.Esacapasepuedeproyectaralaproyecciónortográfica,yaseasobrelamarchaoguardándolaenelCRSortográ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 a 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)
    
respondido por el Jake 23.11.2013 - 22:50
5

Tienes que recortar los datos de tus polígonos a la mitad visible del globo, porque QGIS no lo hace por sí solo.

Escribí un tutorial aquí:

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

EDIT

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

QGIS display world archivos de forma de país centrados en el océano Pacífico con Robinson, Miller Cylindrical u otra proyección

    
respondido por el AndreJ 22.11.2013 - 16:15

Lea otras preguntas en las etiquetas