[OSM-talk] gdalwarp question

[Torsten Mohr]

Hello Jukka,

no, thanks for your help, any hint and discussion is really appreciated.
>
> Sorry, I misunderstood a bit what you were going to do.  It may well be
> that for Mapnik you'll need to reproject raster image first. I do not much
> about Mapnik and while a have been using gdal utilities very much I don't
> believe I have ever needed to create output in Google projection.
> By first look what you have done does make sense. You can make a couple of
> further tests:

I solved the issue now in a different way.  Based on information on wikipedia
about mercaator projection i did the reprojection of the blue marble myself.
It took quite a while to execute, but the reprojection itself worked quite
fine.

Now the borders (shoreline_300) are exactly on the borders of the continents
and islands on the "blue marble".

The relevant part of the script i use is attached below, if anybody is
interested please write me an email (i use an own library to handle the
images themselves, that part needs to be adapted to e.g. PyImage).

Thank you all for providing me help, i hope the script below is of value to
anybody.


Best regards,
Torsten.


#! /usr/bin/python
#coding: latin-1

import ri
from math import *
import gc
import sys

#he = 21600
he = 43200
np = 8
skp_a = 90.0 - (atan(sinh(pi)) * 180.0 / pi)
print skp_a
skp_p = int(he * skp_a / 90.0)
print skp_p

#sys.exit(0)

pat = '/local/vid/earth/strs_%i.tif'

max_phi = atan(sinh(pi))
print "max_phi", max_phi * 180 / pi


def do_img(img, i):
    wi = img.width()
    hi = img.height()
    yi0 = hi / 2

    wo = wi
    ho = he
    yo0 = ho / 2

    print "i", i, "wo", wo, "ho", ho

    out = ri.Ri(wo, ho, 1, 0)

    for x in range(wo):
        for ydo in range(ho / 2):
            yn = float(ydo) / (ho / 2) * pi

            phi = asin(tanh(yn))
            ydi = phi / max_phi * (hi / 2)

            r, g, b, a = img.get(x, int(yi0 + ydi))
            out.set(x, int(yo0 + ydo), r, g, b)

            r, g, b, a = img.get(x, int(yi0 - ydi))
            out.set(x, int(yo0 - ydo), r, g, b)

    out.saveAsPng("/local/vid/earth/ta_%i.png" % (i))
    del out


for i in range(4, np):
    img = ri.Ri(pat % (i), ri.TIFF)
    do_img(img, i)
    del img
    gc.collect()