raster_transformation.go

package main

// export PROJ_LIB=$GOPATH/src/github.com/xeonx/proj4/share/

import (
	"fmt"
	"github.com/xeonx/geom"
	"github.com/xeonx/proj4"
	"image"
	"image/color/palette"
	"image/gif"
	"image/png"
	"log"
	"math"
	"os"
)

const deg2Rad = math.Pi / 180.0
const rad2Deg = 180.0 / math.Pi
const size = 500
const lonStep = 3

type GeoImage struct {
	*image.NRGBA
	X0       float64
	Y0       float64
	PixSizeX float64
	PixSizeY float64
}

func (g *GeoImage) GetLocation(x, y int) (float64, float64) {
	return g.X0 + float64(x)*g.PixSizeX, g.Y0 + float64(y)*g.PixSizeY
}

func (g *GeoImage) GetPixel(x, y float64) (int, int) {
	return int((x - g.X0) / g.PixSizeX), int((y - g.Y0) / g.PixSizeY)
}

func main() {
	projWGS84, err := proj.InitPlus("+init=epsg:4326")
	if err != nil {
		log.Fatal(err)
	}

	infile, err := os.Open("/Users/pablo/Downloads/world.png")
	if err != nil {
		// replace this with real error handling
		panic(err)
	}
	defer infile.Close()

	src, err := png.Decode(infile)
	srcRGBA := src.(*image.NRGBA)
	fmt.Println(srcRGBA.Bounds())
	blueMarble := GeoImage{srcRGBA, -180, 90, 360.0 / 5400.0, -180.0 / 2700.0}

	palette := palette.Plan9
	r := image.Rect(0, 0, size, size)
	var images []*image.Paletted

	for lon := 360.0; lon > 0.0; lon -= lonStep {

		projString := fmt.Sprintf("+proj=ortho +a=6378137.0 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0 +lat_0=%f +lon_0=%f", 0.0, lon)
		fmt.Println(projString)

		projOrtho, err := proj.InitPlus(projString)
		if err != nil {
			log.Fatal(err)
		}

		inv, err := proj.NewTransformation(projOrtho, projWGS84)
		if err != nil {
			log.Fatal(err)
		}

		canvas := GeoImage{image.NewNRGBA(r), -6378127, 6378127, 6378127.0 * 2.0 / float64(r.Dx()), 6378127.0 * -2.0 / float64(r.Dx())}

		globe := make([]geom.Point, r.Dx()*r.Dy())

		for x := 0; x < canvas.Bounds().Dx(); x++ {
			for y := 0; y < canvas.Bounds().Dy(); y++ {
				i, j := canvas.GetLocation(x, y)
				globe[(x*r.Dx())+y] = geom.Point{X: i, Y: j}
			}
		}

		if err := inv.TransformPoints(globe); err != nil {
			log.Fatal(err)
		}

		for x := 0; x < canvas.Bounds().Dx(); x++ {
			for y := 0; y < canvas.Bounds().Dy(); y++ {
				i, j := blueMarble.GetPixel(globe[(x*r.Dx())+y].X*rad2Deg, globe[(x*r.Dx())+y].Y*rad2Deg)
				canvas.Set(x, y, blueMarble.At(i, j))
			}
		}

		paletted := image.NewPaletted(r, palette)
		for i := 0; i < r.Dx(); i++ {
			for j := 0; j < r.Dy(); j++ {
				paletted.Set(i, j, canvas.At(i, j))
			}
		}

		images = append(images, paletted)

	}

	projString := fmt.Sprintf("+proj=ortho +a=6378137.0 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0 +lat_0=%f +lon_0=%f", -28.0, 140.0)

	projOrtho, err := proj.InitPlus(projString)
	if err != nil {
		log.Fatal(err)
	}

	inv, err := proj.NewTransformation(projOrtho, projWGS84)
	if err != nil {
		log.Fatal(err)
	}

	//Save the animated GIF file
	imgFile, err := os.Create("rotating_earth.gif")
	if err != nil {
		log.Fatal(err)
	}

	delays := make([]int, len(images))
	for i := range delays {
		delays[i] = 1
	}

	g := gif.GIF{
		Image: images,
		Delay: delays,
	}

	err = gif.EncodeAll(imgFile, &g)
	if err != nil {
		log.Fatal(err)
	}

	canvas := GeoImage{image.NewNRGBA(r), -6378127, 6378127, 6378127.0 * 2.0 / float64(r.Dx()), 6378127.0 * -2.0 / float64(r.Dx())}

	globe := make([]geom.Point, r.Dx()*r.Dy())

	for x := 0; x < canvas.Bounds().Dx(); x++ {
		for y := 0; y < canvas.Bounds().Dy(); y++ {
			i, j := canvas.GetLocation(x, y)
			globe[(x*r.Dx())+y] = geom.Point{X: i, Y: j}
		}
	}

	if err := inv.TransformPoints(globe); err != nil {
		log.Fatal(err)
	}

	for x := 0; x < canvas.Bounds().Dx(); x++ {
		for y := 0; y < canvas.Bounds().Dy(); y++ {
			i, j := blueMarble.GetPixel(globe[(x*r.Dy())+y].X*rad2Deg, globe[(x*r.Dy())+y].Y*rad2Deg)
			canvas.Set(x, y, blueMarble.At(i, j))
		}
	}

	out, err := os.Create("earth.png")
	if err != nil {
		fmt.Println(err)
		os.Exit(1)
	}
	png.Encode(out, canvas)

	fmt.Println("DONE!")
}