Unable to use GPS coordinates for GEOTiff map

Question

I have been using a map to calculate some distances at sea. The map had a bit low resolution so I got another one with better resolution (both GEOTiff), but with the amounts of distances I need to calculate the algorithm quickly grew to require a run time of 6 months to complete. I want to try to make cropped versions of the map with the point I am calculating distances for in the center, and a box with sides of 200km around it. I found some example code using rasterio to do this:

MIN_LAT = 57.75
MAX_LAT = 71.35
MIN_LON = 6
MAX_LON = 31.45
r = 200
c = (22.3542, 70.14145)

dlon, dlat = 150*0.015060, 100*0.008983


def check_val(val, latorlon):
    if latorlon == "lat":
        if val < MIN_LAT:
            return MIN_LAT
        elif val > MAX_LAT:
            return MAX_LAT
    elif latorlon == "lon":
        if val < MIN_LON:
            return MIN_LON
        elif val > MAX_LON:
            return MAX_LON

    return val

points = [(check_val(c[0] - dlon, "lon"), check_val(c[1] - dlat, "lat")),
           (check_val(c[0] + dlon, "lon"), check_val(c[1] - dlat, "lat")),
           (check_val(c[0] + dlon, "lon"), check_val(c[1] + dlat, "lat")),
           (check_val(c[0] - dlon, "lon"), check_val(c[1] + dlat, "lat"))]

geoms = [{'type': 'Polygon', 'coordinates': [points]}]

with rasterio.open('map_100x100_original.tif') as src:
    out_image, out_transform = mask(src, geoms, crop=True)
out_meta = src.meta.copy()

out_meta.update(
    {
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    }
)

with rasterio.open("map_100x100_cropped.tif", "w", **out_meta) as dest:
    dest.write(out_image)

The code works as intended with the old map, however, with the new map I get an error message from rasterio saying: "WindowError: windows do not intersect" followed by "ValueError: Input shapes do not overlap raster".

Using some code I found here I got the following information on the two maps:

old map

None
WGS 84

new map

WGS 84 / UTM zone 33N
WGS 84

Using some other code I found here I got the corners of the two maps

old map

(-180.0, 90.00000000000001) (180.00000000007202, -90.000000000036)

new map

(1121948.79, 6426051.97) (-99551.21, 7962751.97)

Using QGIS I got some coordinates by hovering on the map which I could plug in to the code and it would work. For example, the GPS coordinates below

(59.48129, 5.89213)

Would be something like

(-60045, 6609466)

But I need to be able to use GPS coordinates, so is there a way I could fix the new map?

score 1 · Accepted Answer · answered Oct 03 '18 at 12:43

I was able to solve my problem using this code:

raster = gdal.Open("url/to/map.tif")
srs = osr.SpatialReference()
srs.ImportFromWkt(raster.GetProjection())
srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srsLatLong, srs)
x, y, height = ct.TransformPoint(lon, lat)

and then using the x and y in stead of the latitude and longitude coordinates

score 0 · Answer 2 · answered Oct 03 '18 at 12:42

I used your "new map points" (WGS 84 / UTM zone 33N; EPSG: 32633):

(1121948.79, 6426051.97) (-99551.21, 7962751.97)

for generating a bounding box for this area. Afterward, I got a map by using an own plugin to clip by mask referred area. Coordinates are in WGS 84 / UTM zone 33N. It looks as:

So, for producing a map with coordinates of your "old map" you need to reproject to WGS 84, EPSG: 4326.

If you want to use it in python you can adapt following command:

gdalwarp -overwrite -s_srs EPSG:32633 -of GTiff /home/zeito/Desktop/map_32633_clipped.tif /home/zeito/Desktop/map_4326_clipped.tif

After running above command, map in EPSG 4326 coordinates looks as:

By using a UTM_Zone_Boundaries shapefile it can be observed that WGS 84 / UTM zone 33N is colored in yellow.

You can clip your map for this overlapped area.

Thank you for your answer. I was able to solve my problem using coordinate transformations in my code, so I don't have to change my map after all :) — niknoe, Oct 03 '18 at 12:45

Unable to use GPS coordinates for GEOTiff map

2 Answers2