1

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?

1

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

0

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:

enter image description here

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:

enter image description here

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

enter image description here

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 :) – Niklas Noem Oct 3 '18 at 12:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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