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.
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?
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)
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:
