EOxServer / eoxserver

EOxServer is a Python application and framework for presenting Earth Observation (EO) data and metadata.
https://eoxserver.org
Other
40 stars 19 forks source link

Registration for images with transposed GeoTransform currently works but subsequent requests fail #439

Open Schpidi opened 4 years ago

Schpidi commented 4 years ago

Because we're always using elements 1 and 5 of the GeoTransform for the offset vector which are 0 in the case of transposed images the extent is 0 which yields an error.

Ideally we'd support serving images like this but as long as MapServer can't handle them we should fail during ingestion with a proper failure message.

lubojr commented 4 years ago

Or fix the image geotransform during registration:

def correct_geo_transform(src_dst):
    ulx, xres, xskew, uly, yskew, yres = src_dst.GetGeoTransform()
    # test geotransform if necessary to shift
    if xres == 0.0 and yres == 0.0:
        # malformed image, compute xres and yres switched in geotransform
        lrx = ulx + (src_dst.RasterXSize * xskew)
        lry = uly + (src_dst.RasterYSize * yskew)
        # [ulx, lrx, lry, uly] - bounds = lon_min, lon_max, lat_min, lat_max
        fp = [[0, src_dst.RasterXSize, src_dst.RasterXSize, 0], [0, 0, src_dst.RasterYSize, src_dst.RasterYSize]]
        tp = [[ulx, lrx, lrx, ulx], [lry, lry, uly, uly]]
        pix = list(zip(fp[0], fp[1]))
        coor = list(zip(tp[0], tp[1]))
        # compute the gdal.GCP parameters
        gcps = []
        for index, txt in enumerate(pix):
            gcps.append(gdal.GCP())
            gcps[index].GCPPixel = pix[index][0]
            gcps[index].GCPLine = src_dst.RasterYSize - int(pix[index][1])
            gcps[index].GCPX = coor[index][0]
            gcps[index].GCPY = coor[index][1]
        # get correct geotransform from gcps
        geotransform_new = gdal.GCPsToGeoTransform(gcps)
        # overwrite geotransform with new
        src_dst.SetGeoTransform(geotransform_new)
    return src_dst