5

Image corner coordinates are known:

  • Top left: (70.20315, 20.92749)
  • Top right: (70.20315, 31.45240)
  • Bottom left: (59.70571, 20.92749)
  • Bottom right: (59.70571, 31.45240)

(Basically map over Finland and a little bit of Sweden/Norway/Russia)

I am not sure what is current projection name; +1 in latitude and longitude is constant distance in image (which is not the case with Mercator projection).

I couldn't find a way to change image projection to Google maps overlay image, preferably using Python. Python isn't mandatory, but command line interface is.

It's possible to do that with pyproj and PIL, but I wasn't able to figure out how to do that. I also found GDAL, and more specifically gdalwarp, but couldn't find options for changing projection in a correct way.

Example image

4 Answers 4

3

One solution is to use gdalwarp, as -s_srs and -t_srs options you can use the definitions or EPSG codes found in this site.

Since you don't know the source srs you can try some options and check the results with another reliable data. But is more recommended that you trace the source of this image and try to find it's srs definition.

If the image is a geotiff, it may already contain the srs definition. You can check it with gdalinfo. GDAL.

2
  • Well, source of that image is another python code, which calculates heatmap. Changing that code is much more problematic. It basically just calculates constant grid over the area, without any projections. So can't find any help from there. And no srs definitions included (obviously).
    – Olli
    Commented Jun 1, 2011 at 7:30
  • @Olli: The data used as input into the heatmap script has to be in some projection. Do you know the source? What exactly are you interpolating here? It's obvious that the heatmap script fails on the image border areas.
    – underdark
    Commented Jun 1, 2011 at 8:46
2

Because your corner coordinates form a square and your image is also square, I would say that it was in a "Plate Carree"/equidistant cylindrical projection, which corresponds to EPSG:32663. So a gdalwarp call should read something like:

gdalwarp -s_srs EPSG:32663 -t_srs EPSG:4326 -r cubic srcimage.tif destimage.tif

Choose your own -t_srs parameter if WGS84 isn't what you want, in which case you'll probably need a -dstnodata option, and you may need to play aorund with different -r parameters.

0

Your description of the projection suggests it is Equirectangular as merseyviking said

0

Thanks for answers. However, extensive trial-and-error showed it worked fine without specifying source projection. Steps to reproduce:

Adding control points to TIF (of course, first convert mapimage.png mapimage.tif):

gdal_translate mapimage.tif mapimage_georeffed.tif -gcp 2300 2300 31.45240 59.70571 -gcp 0 0 20.92749 70.20315 -gcp 2300 0 31.45240 70.20315 -gcp 0 2300 20.92749 59.70571 -a_srs +proj=longlat

Then, GDAL python bindings include gdal2tiles.py, which autogenerates tiles on right projection for maps, and also sample map pages.

For the record, I'm not absolutely positive projection is exactly correct, but over that area it seems to be at least good enough.

2
  • gdal_translate will only insert GCPs (ground control points) into your image, it won't actually do any warping. I've used your method in the past with earth resistivity data that was rotated away from north, but I didn't want any interpolation in the data until it came to display it. It seems that proj.4 assumes WGS84 ellipsoid and datum, but I would be inclined to add +ellps=WGS84 +datum=WGS84 +nodefs to the end to make sure. Commented Jun 2, 2011 at 9:30
  • Yes, gdal_translate do not warp image, but gdal2tiles.py does. In the same time it handles tiling, which is added bonus for goal (having image as working overlay in interactive map).
    – Olli
    Commented Jun 2, 2011 at 9:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.