You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 

215 lines
8.0 KiB

# Builds simple Tissot Indicatrices from the cell centres of a DGG.
#
# This file needs the Homolosine projection counter domain. It can be obtained
# from Zenodo: https://zenodo.org/record/1841337
#
# Author: Luís de Sousa luis (dot) de (dot) sousa (at) protonmail (dot) ch
# Date: 24-11-2018
###############################################################################
import math
import numpy as np
import matplotlib.pyplot as plt
from geographiclib.geodesic import Geodesic
from osgeo import ogr
from osgeo import osr
import shapely.wkb
from shapely.geometry import LineString
import csv
DEBUG = False
radius = 1000 # metres
azimuths = [0, 90, 180, 270] # In GeographicLib equivalent to [N,E,S,W]
#centresFile = "data/centresLand.gpkg"
#pathOutput = "data/land/"
centresFile = "data/centres.gpkg"
pathOutput = "data/globe/"
counterDomainFile = "data/CounterDomain.geojson"
counterDomain = None
geographic = osr.SpatialReference()
homolosine = osr.SpatialReference()
hammer = osr.SpatialReference()
sinusoidal = osr.SpatialReference()
eckertIV = osr.SpatialReference()
mollweide = osr.SpatialReference()
geographic.ImportFromProj4("+proj=longlat +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs")
homolosine.ImportFromProj4("+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs")
hammer.ImportFromProj4("+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs +wktext")
sinusoidal.ImportFromProj4("+proj=sinu +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs")
eckertIV.ImportFromProj4("+proj=eck4 +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs")
mollweide.ImportFromProj4("+proj=moll +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs")
transformHomolosine = osr.CoordinateTransformation(geographic, homolosine)
transformHammer = osr.CoordinateTransformation(geographic, hammer)
transformSinusoidal = osr.CoordinateTransformation(geographic, sinusoidal)
transformEckertVI = osr.CoordinateTransformation(geographic, eckertIV)
transformMollweide = osr.CoordinateTransformation(geographic, mollweide)
def getIndicatricesCentres():
print("Getting ready to load centres layer")
driver = ogr.GetDriverByName("GPKG")
dataSource = driver.Open(centresFile, 0)
return dataSource
def loadHomolosineCounterDomain():
global counterDomain
driver = ogr.GetDriverByName("GeoJSON")
dataSource = driver.Open(counterDomainFile, 0)
layer = dataSource.GetLayer()
feature = layer.GetFeature(0)
counterDomain = shapely.wkb.loads(feature.GetGeometryRef().ExportToWkb())
def computeAngle(centre, point, az):
inverse = 1
diffX = point.GetPoint(0)[0] - centre.GetPoint(0)[0]
diffY = point.GetPoint(0)[1] - centre.GetPoint(0)[1]
if (az == azimuths[0] or az == azimuths[2]):
adjacent = diffY
opposed = diffX
else:
adjacent = diffX
opposed = diffY
inverse = -1
if adjacent == 0:
return 0
else:
return inverse * math.degrees(math.atan(opposed / adjacent))
def computeDistance(centre, point):
line = LineString([(centre.GetPoint(0)[0], centre.GetPoint(0)[1]),
(point.GetPoint(0)[0], point.GetPoint(0)[1])])
return line.length - radius
def computeDistanceHomolosine(centre, point):
line = LineString([(centre.GetPoint(0)[0], centre.GetPoint(0)[1]),
(point.GetPoint(0)[0], point.GetPoint(0)[1])])
if(line.length > 1500):
intersection = line.intersection(counterDomain)
return abs(intersection.length - radius)
else:
return abs(line.length - radius)
def computeDistortions(layer, geoTransform, distanceFunction):
errors = {"angle":[], "distance":[]}
errorsIndicatrix = {"long":[], "lat":[], "angle":[], "distance":[]}
for feature in layer:
geom = feature.GetGeometryRef()
centre = ogr.Geometry(ogr.wkbPoint)
if(DEBUG):
print("The point recieved: ", str(geom.GetX()), " ", str(geom.GetY()))
centre.AddPoint(geom.GetX(), geom.GetY())
centre.Transform(geoTransform)
if(DEBUG):
print("The point transformed: ", str(centre.GetX()), " ", str(centre.GetY()))
absErrorAngle = 0
absErrorDistance = 0
for az in azimuths:
# Latitude comes first !
geodesic = Geodesic.WGS84.Direct(geom.GetY(), geom.GetX(), az, radius)
vertex = ogr.Geometry(ogr.wkbPoint)
vertex.AddPoint(geodesic["lon2"], geodesic["lat2"])
vertex.Transform(geoTransform)
if(DEBUG):
print("Vertex ", str(az), ": ", str(vertex.GetX()), " ", str(vertex.GetY()))
errorAngle = computeAngle(centre, vertex, az)
errorDistance = distanceFunction(centre, vertex)
errors["angle"].append(errorAngle)
errors["distance"].append(errorDistance)
absErrorAngle += abs(errorAngle)
absErrorDistance += abs(errorDistance)
errorsIndicatrix["long"].append(geom.GetX())
errorsIndicatrix["lat"].append(geom.GetY())
errorsIndicatrix["angle"].append(absErrorAngle / 4)
errorsIndicatrix["distance"].append(absErrorDistance / 4)
if(DEBUG):
input("Continue?")
return errors,errorsIndicatrix
def showSummaryStatistics(results, name):
print("\n" + name)
print("==============")
print("---- Angles ----")
print("Mean: " + str(np.mean(results["angle"])))
print("Variance: " + str(np.var(results["angle"])))
print("Standard deviation: " + str(np.std(results["angle"])))
print("RMSE: " + str(np.sqrt(np.mean(np.power(results["angle"], 2)))))
print("MAE: " + str(np.mean(np.abs(results["angle"]))))
#print("R-squared: " + )
print("---- Distances ----")
print("Mean: " + str(np.mean(results["distance"])))
print("Variance: " + str(np.var(results["distance"])))
print("Standard deviation: " + str(np.std(results["distance"])))
print("RMSE: " + str(np.sqrt(np.mean(np.power(results["distance"], 2)))))
print("MAE: " + str(np.mean(np.abs(results["distance"]))))
print("Number of measurements: " + str(len(results["distance"])))
def writeOutput(errors, errorsIndicatrix, name):
with open(pathOutput + name + "_mapping.csv", 'w') as myfile:
wr = csv.writer(myfile, quoting=csv.QUOTE_ALL)
wr.writerow(('Longitude', 'Latitude', 'Angle', 'Distance'))
for i in range(0, len(errorsIndicatrix["angle"])):
wr.writerow((errorsIndicatrix["long"][i], errorsIndicatrix["lat"][i],
errorsIndicatrix["angle"][i], errorsIndicatrix["distance"][i]))
with open(pathOutput + name + ".csv", 'w') as myfile:
wr = csv.writer(myfile, quoting=csv.QUOTE_ALL)
wr.writerow(('Angle', 'Distance'))
for i in range(0, len(errors["angle"])):
wr.writerow((errors["angle"][i], errors["distance"][i]))
def main():
dataSource = getIndicatricesCentres()
layer = dataSource.GetLayer(0)
print("Number of Inidicatrices to process: ", layer.GetFeatureCount())
loadHomolosineCounterDomain()
errors, errorsIndicatrix = computeDistortions(layer, transformHomolosine, computeDistanceHomolosine)
writeOutput(errors, errorsIndicatrix, "Homolosine")
showSummaryStatistics(errors, "Homolosine")
errors, errorsIndicatrix = computeDistortions(layer, transformSinusoidal, computeDistance)
writeOutput(errors, errorsIndicatrix, "Sinusoidal")
showSummaryStatistics(errors, "Sinusoidal")
errors, errorsIndicatrix = computeDistortions(layer, transformEckertVI, computeDistance)
writeOutput(errors, errorsIndicatrix, "EckertIV")
showSummaryStatistics(errors, "EckertIV")
errors, errorsIndicatrix = computeDistortions(layer, transformMollweide, computeDistance)
writeOutput(errors, errorsIndicatrix, "Mollweide")
showSummaryStatistics(errors, "Mollweide")
errors, errorsIndicatrix = computeDistortions(layer, transformHammer, computeDistance)
writeOutput(errors, errorsIndicatrix, "Hammer")
showSummaryStatistics(errors, "Hammer")
if __name__ == '__main__': main()