# Unable to use GPS coordinates for GEOTiff map

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?

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

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 :) Oct 3, 2018 at 12:45