I`m trying to create a sphere-like view using qgis and the "world from space"-projection http://spatialreference.org/ref/sr-org/6980/ (essentialy an ortho-projection). ArcGIS wraps the shapes correctly but QGIS (2.01) produces nasty artifacts.

enter image description here

I have to produce the globes on a regular basis with different angles so does anyone out there have an idea how to fix this problem?

  • 1
    related QGIS bug report: hub.qgis.org/issues/2703 – naught101 Jan 6 '15 at 7:05
  • Is it too big a technical problem to have an orthographic projection pre-loaded, that can be re-centred to any view? – user66708 Feb 6 '16 at 13:42
  • This does not answer the question. Please take the tour to learn how to ask a focussed question. – John Powell aka Barça Feb 6 '16 at 13:49
up vote 22 down vote
+200

As Andre said, for this to work, you'll need to crop your layer before projecting it. Andre describes a manual method, which works well for a lot of cases: Project your shapefile to an azimuthal equidistant projection with the same parameters as the orthographic projection, create a clipping circle that covers the hemisphere that will be visible in the orthographic projection, and clip the shapefile with that. However, that method requires a fair bit of manual effort, and doesn't work for all projection parameters, since projecting to an azimuthal equidistant projection can lead to similar problems as projecting to an orthographic projection.

Here's a script (now also available as the Clip to Hemisphere QGIS plugin) that takes a slightly different approach: A clipping layer is created in the coordinate reference system of the original shapefile by projecting a circle from the orthographic to the source CRS, but additionally making sure to cover the whole visible hemisphere, including the visible pole.

This is what the clipping layer looks like for an orthographic projection centered on 30°N, 110°E:

The script then clips the currently selected layer with the clipping layer, and adds the resulting layer to the project. That layer can then be projected to the orthographic projection, either on the fly or by saving it in the orthographic CRS:

Here's the script. Make sure to save it in your Python path, for example as 'cliportho.py'. Then you can import it in the QGIS Python console using import cliportho. To clip a layer, call 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)
  • Looks very promising - I will definitely try this out and be happy to provide feedback. I`m a little into arcpy programming but have not started with qgis programming - but I will try to understand what you are doing ;-) A plugin (maybe working batch for several layers) would be so helpful! – user1523709 Nov 24 '13 at 9:34
  • 2
    FWIW this works exceptionally well. Thanks! – Hamish Nov 20 '14 at 19:49
  • 1
    FYI, This script no longer works in QGIS 2.16, due to the removal of the "fTools" package. – Spike Williams Aug 21 '16 at 21:30
  • 2
    @SpikeWilliams: I've updated the script to remove the dependence on fTools. – Jake Aug 22 '16 at 19:37
  • 1
    Works great! Thanks! – Spike Williams Aug 24 '16 at 1:22

You have to crop your polygon data to the visible half of the globe, because QGIS does not do that by itself.

I wrote a tutorial here:

Where did the polygons go after projecting a map in QGIS?


EDIT

The picture you show is actually not an ortho projection, as it shows the whole world, and not only the visible half as seen from outer space. For world maps, the cutting is a bit easier, as described here:

QGIS display world country shape files centered on pacific ocean using Robinson, Miller Cylindrical or other projection

  • Thanks Andre, that was quite helpful to understand the problem - but since I have to create such globes on a near-daily basis and with changing perspectives it requires a lot of manual work. Do you know of any plugins-etc. to automate your solution? – user1523709 Nov 22 '13 at 15:35
  • Once you created a clipping circle, the rest can be done using GDAL at the command line level using a batch script. – AndreJ Nov 24 '13 at 6:56

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.