library(rgdal)
library(raster)	# for crop()
library(maptools) # for data(wrld_simpl)

# http://proj4.org/geodesic.html
wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

world_rectangle = matrix(ncol=2, byrow=T, data=c(-180, 90, 180, 90, 180, -90, -180, -90, -180, 90))
world_rectangle = SpatialPolygons(list(Polygons(list(Polygon(world_rectangle)), 1)), proj4string = wgs84)

graticule = gridlines(world_rectangle, easts = seq(-180, 180, 10), norths = seq(-90, 90, 10))

data(wrld_simpl)
countries = spTransform(wrld_simpl, wgs84)


# Plotting Function

mmsp_plot = function(countries, graticule, projection, outfile, height=400, width=600) {
	if (class(projection) == "CRS") {
		countries = spTransform(countries, projection)
		graticule = spTransform(graticule, projection) }

	print(graticule)

	png(filename = outfile, width=width, height=height)

	par(mar=c(0,0,0,0))

	# Plot the graticule first so edges are not cut off
	plot(graticule, col="white")
	plot(countries, add=T, border="#c0c0c0", col="#c0c0c0")
	plot(graticule, add=T, col="#808080")

	dev.off()
}

# Convenience function for cropping Spatial* to a rectangle

mmsp_crop = function(objects, left, top, right, bottom) {
	wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
	objects = spTransform(objects, wgs84)
	cropbox = matrix(ncol=2, byrow=T, data=c(left, top, right, top, right, bottom, left, bottom, left, top))
	cropbox = SpatialPolygons(list(Polygons(list(Polygon(cropbox)), 1)), proj4string = wgs84)
	return(crop(objects, cropbox))
}

mmsp_apply = function(sp, fun) {
	if (class(sp) == "SpatialLines") {
		for (z1 in 1:length(sp@lines)) {
			for (z2 in 1:length(sp@lines[[z1]]@Lines)) {
				for (z3 in 1:nrow(sp@lines[[z1]]@Lines[[z2]]@coords)) {
					sp@lines[[z1]]@Lines[[z2]]@coords[z3,] = 
						fun(sp@lines[[z1]]@Lines[[z2]]@coords[z3,])
					}
				}
			}
		return(SpatialLines(sp@lines, proj4string=sp@proj4string))
	}

	if (class(sp) == "SpatialPolygonsDataFrame") {
		for (z1 in 1:length(sp@polygons)) {
			for (z2 in 1:length(sp@polygons[[z1]]@Polygons)) {
				for (z3 in 1:nrow(sp@polygons[[z1]]@Polygons[[z2]]@coords)) {
					sp@polygons[[z1]]@Polygons[[z2]]@coords[z3,] = 
						fun(sp@polygons[[z1]]@Polygons[[z2]]@coords[z3,])
					}
				}
			}
		return(SpatialPolygonsDataFrame(SpatialPolygons(sp@polygons, proj4string=sp@proj4string), sp@data, match.ID = F))
	}
}


# Scales to stretch, compress, or invert

mmsp_scale = function(sp, scale_x = 1, scale_y = 1) {
	scalefun = function(z) return(matrix(data=c(z[1] * scale_x, z[2] * scale_y), ncol=2))
	return(mmsp_apply(sp, scalefun))
}


# Recenters Spatial* to an arbitrary center, creates WGS 84 false latitude and longitude

mmsp_recenter = function(objects, center_x = 0, center_y = 0) {
	if (center_x < 0) {
		xcut = center_x + 180 
		o1 = mmsp_crop(objects, -180, 90, xcut, -90)
		o1 = shift(o1, x = -center_x)
		o2 = mmsp_crop(objects, xcut, 90, 180, -90)
		o2 = shift(o2, x = -360 - center_x)
		objects = rbind(o1, o2, makeUniqueIDs = T)
	} else if (center_x > 0) {
		xcut = center_x - 180 
		o1 = mmsp_crop(objects, xcut, 90, 180, -90)
		o1 = shift(o1, x = -center_x)
		o2 = mmsp_crop(objects, -180, 90, xcut, -90)
		o2 = shift(o2, x = 360 - center_x)
		objects = rbind(o1, o2, makeUniqueIDs = T)
	}

	if (center_y > 0) {
		ycut = center_y - 90
		o1 = mmsp_crop(objects, -180, 90, 180, ycut)
		o1 = shift(o1, y = -center_y)
		o2 = mmsp_crop(objects, -180, 90, 180, 90 - center_y)
		o2 = mmsp_scale(o2, -1, -1)
		o2 = shift(o2, y = 180 - center_y)
		objects = rbind(o1, o2, makeUniqueIDs = T)
	} else if (center_y < 0) {
		ycut = center_y + 90
		o1 = mmsp_crop(objects, -180, ycut, 180, -90)
		o1 = shift(o1, y = -center_y)
		o2 = mmsp_crop(objects, -180, -90 - center_y, 180, 90)
		o2 = mmsp_scale(o2, -1, -1)
		o2 = shift(o2, y = -180 - center_y)
		objects = rbind(o1, o2, makeUniqueIDs = T)
	}

	return(objects)
}

# Transform to side-by-side hemispheres
mmsp_hemispheres = function(objects, central_meridian1, central_meridian2, crs) {
	wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
	objects = spTransform(objects, wgs84)

	robjects1 = mmsp_recenter(objects, central_meridian1)
	robjects1 = mmsp_crop(robjects1, -90, 90, 90, -90)
	robjects1 = spTransform(robjects1, crs)

	robjects2 = mmsp_recenter(objects, central_meridian2)
	robjects2 = mmsp_crop(robjects2, -90, 90, 90, -90)
	robjects2 = spTransform(robjects2, crs)

	robjects2 = shift(robjects2, x = robjects1@bbox[1,2] - robjects2@bbox[1,1])

	return(rbind(robjects1, robjects2))
}

# Azimuthal Equal-Area
# http://spatialreference.org/ref/sr-org/7250/
azimuthal_equal_area_equatorial = CRS("+proj=laea +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
rcountries = mmsp_crop(countries, -160, 90, 160, -90)
rgraticule = mmsp_crop(graticule, -160, 90, 160, -90)
mmsp_plot(rcountries, rgraticule, azimuthal_equal_area_equatorial, "azimuthal-equal-area-equatorial.png")

azimuthal_equal_area_polar = CRS("+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
rcountries = mmsp_crop(countries, -180, 80, 180, -60)
rgraticule = mmsp_crop(graticule, -180, 80, 180, -80)
mmsp_plot(rcountries, rgraticule, azimuthal_equal_area_polar, "azimuthal-equal-area-polar.png")


# Azimuthal Equidistant, Polar
azimuthal_equidistant_polar = CRS("+proj=aeqd +lat_0=90 +lon_0=-70")
rcountries = mmsp_crop(countries, -180, 90, 180, -70)
mmsp_plot(rcountries, graticule, azimuthal_equidistant_polar, "azimuthal-equidistant-polar.png")

# Azimuthal Equidistant, US-Centered
azimuthal_equidistant_us = CRS("+proj=aeqd +lat_0=40 +lon_0=-100")
mmsp_plot(countries, graticule, azimuthal_equidistant_us, "azimuthal-equidistant-us.png")


# EPSG:800400 Aitoff = stretched azimuthal equidistant
# Aitoff Projection with central Meridian 90°W. Modified azimuthal. Neither conformal nor equal area

aitoff = CRS("+proj=aitoff")
mmsp_plot(countries, graticule, aitoff, "aitoff.png")

# Hammer projection = stretched azimuthal equal-area
mmsp_plot(countries, graticule, hammer, "hammer.png")
stop()




# Briesemeister = Oblique Hammer variant
# https://download.osgeo.org/proj/proj.4.3.I2.pdf
briesemeister = CRS("+proj=ob_tran +o_proj=hammer +o_lat_p=45 +o_lon_p=0 +lon_0=10 +M=0.93541")
rgraticule = mmsp_crop(graticule, -170, 80, 175, -80)
rcountries = mmsp_crop(countries, -180, 90, 180, -70)
mmsp_plot(rcountries, rgraticule, briesemeister, "briesemeister.png")


# Globular - Apian

apian = CRS("+proj=apian")
rcountries = mmsp_hemispheres(countries, -110, 50, apian)
rgraticule = mmsp_hemispheres(graticule, -110, 50, apian)
mmsp_plot(rcountries, rgraticule, NULL, "globular-apian.png")

rcountries = mmsp_hemispheres(countries, -170, -20, apian)
rgraticule = mmsp_hemispheres(graticule, -170, -20, apian)
mmsp_plot(rcountries, rgraticule, NULL, "globular-apian-oceans.png")


# Globular - Bacon

bacon = CRS("+proj=bacon")
rcountries = mmsp_hemispheres(countries, -110, 50, bacon)
rgraticule = mmsp_hemispheres(graticule, -110, 50, bacon)
mmsp_plot(rcountries, rgraticule, NULL, "globular-bacon.png")

# Globular - Nicolosi

nicolosi = CRS("+proj=nicol")
rcountries = mmsp_hemispheres(countries, -110, 50, nicolosi)
rgraticule = mmsp_hemispheres(graticule, -110, 50, nicolosi)
mmsp_plot(rcountries, rgraticule, NULL, "globular-nicolosi.png")

# Goode Homolosine
# Goode's homolosine map projection is designed to minimize distortion 
# for the entire world. It is an interrupted pseudocylindrical equal-area projection. 
# John Paul Goode developed the projection in 1925.

goode_mask = matrix(ncol=2, byrow=T, data=c(
-180, 90,
-40, 90, -40, 0, -30, 0, -30, 90,
180, 90, 180, -90,
90, -90, 90, 0, 80, 0, 80, -90,
-10, -90, -10, 0, -20, 0, -20, -90,
-90, -90, -90, 0, -100, 0, -100, -90,
-180, -90, -180, 90))
goode_mask = SpatialPolygons(list(Polygons(list(Polygon(goode_mask)), 1)), proj4string = wgs84)

rcountries = crop(countries, goode_mask)
rgraticule = crop(graticule, goode_mask)

goode_homolocine = CRS("+proj=igh")
mmsp_plot(rcountries, rgraticule, goode_homolocine, "goode_homolosine.png")

# Lagrange
# World at War
# This is the 1902 maps projection?!

lagrange = CRS("+proj=lagrng")
rcountries = mmsp_recenter(countries, -90)
rcountries = mmsp_crop(rcountries, -180, 90, 180, -80)
rgraticule = mmsp_crop(graticule, -180, 90, 180, -90)
mmsp_plot(rcountries, rgraticule, lagrange, "lagrange.png")


# EPSG: 3395 World Mercator
world_mercator = CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
rcountries = mmsp_crop(countries, -180, 80, 180, -80)
rgraticule = mmsp_crop(graticule, -180, 80, 180, -80)
mmsp_plot(rcountries, rgraticule, world_mercator, "mercator.png")

# Mercator: MacArthur's Universal Corrective

world_mercator = CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
rcountries = mmsp_crop(countries, -180, 80, 180, -80)
rgraticule = mmsp_crop(graticule, -180, 80, 180, -80)
rcountries = mmsp_scale(rcountries, 1, -1)

mmsp_plot(rcountries, rgraticule, world_mercator, "mcarthurs-corrective.png")


# Nordic = Oblique Hammer variant
nordic = CRS("+proj=ob_tran +o_proj=hammer +o_lat_p=45 +o_lon_p=0 +lon_0=0 +M=1")
rcountries = mmsp_crop(countries, -180, 90, 180, -70)
mmsp_plot(rcountries, graticule, nordic, "nordic.png")


# Erwin Raisz's Armadillo projection (1943)
# https://www.wired.com/2014/01/projection-raisz-armadillo/

mmsp_raisz_armadillo = function(objects) {
	# https://en.wikipedia.org/wiki/Armadillo_projection
	radius = 6353 # km
	armadillo_transform = function(coord) {
			long = coord[1] * pi / 180
			lat = coord[2] * pi / 180
			twenty = 28 * pi / 180
			x = radius * (1 + cos(lat)) * sin(long / 2)
			y = radius * (((1 + sin(twenty) - cos(twenty)) / 2) + (sin(lat) * cos(twenty)) -
				((1 + cos(lat)) * sin(twenty) * cos(long / 2)))
			return(c(x, y))
		}

	objects = mmsp_apply(objects, armadillo_transform)
	objects.proj4string = CRS("")
	return(objects)
}

rcountries = mmsp_crop(countries, -180, 90, 180, -70)
rgraticule = mmsp_crop(graticule, -180, 90, 180, -70)

rcountries = mmsp_raisz_armadillo(rcountries)
rgraticule = mmsp_raisz_armadillo(rgraticule)

mmsp_plot(rcountries, rgraticule, NULL, "raisz-armadillo.png")
	

# Winkel One
# Elementi Climatologici
winkel1 = CRS("+proj=winkl")
mmsp_plot(countries, graticule, winkel1, "winkel-one.png")

# SR-ORG:7291 Winkel Triple
winkel = CRS("+proj=wintri")
mmsp_plot(countries, graticule, winkel, "winkel-tripel.png")



# EPSG: 54030 World Robinson
robinson = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
mmsp_plot(countries, graticule, robinson, "robinson.png")

# Equirectangular (Unprojected WGS 84)
mmsp_plot(countries, graticule, wgs84, "equirectantular.png")


# EPSG: 3857 Spherical Mercator
spherical_mercator = CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 
	+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs")

mmsp_plot(countries, graticule, spherical_mercator, "spherical_mercator.png")


# Bartholemew’s Regional Projection (Interrupted)
# Fuller's Dymaxion
# Gingery Projection
# Sinu-Mollweide?
# Sinusoidal Equal Area
# Winkel II?