TransWikia.com

gdalwrap incorrectly re-project wide angle

Geographic Information Systems Asked by Andrey Rodionov on December 4, 2020

I’m trying to re-project Meteor-M 2 satellite image with gdalwrap. I have GCP points which I put into .vrt file:

<VRTDataset rasterXSize="1568" rasterYSize="480">
  <GCPList Projection="EPSG:4326">
    <GCP Id="0" Pixel="1568.5" Line="0.5" X="62.20476647354344" Y="42.09684997440528" Z="0.0" />
    ... lots of them
    </GCPList>
    <VRTRasterBand dataType="Byte" band="1">
      <SimpleSource>
        <SourceFilename relativeToVRT="1">output.png</SourceFilename>
        <SourceBand>1</SourceBand>
      </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

Then I executed gdalwrap with the following arguments:

gdalwarp -tps -overwrite -t_srs epsg:3857 -of GTIFF output.vrt output10.tif

Then I cut this geotiff into tiles:

gdal2tiles.py –profile=mercator -z 3-8 output10.tif tiles

So far so good. But when I draw both GCP and tiles onto open street map I’m getting:

meteor-m2

Re-projected image and GCPs differ greatly.

gdalinfo:

Driver: GTiff/GeoTIFF
Files: output10.tif
Size is 1656, 583
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (2926801.391389018390328,6009567.154649123549461)
Pixel Size = (2411.811113223014218,-2411.811113223014218)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2926801.391, 6009567.155) ( 26d17'30.86"E, 47d24'42.83"N)
Lower Left  ( 2926801.391, 4603481.276) ( 26d17'30.86"E, 38d10'12.31"N)
Upper Right ( 6920760.595, 6009567.155) ( 62d10'12.90"E, 47d24'42.83"N)
Lower Right ( 6920760.595, 4603481.276) ( 62d10'12.90"E, 38d10'12.31"N)
Center      ( 4923780.993, 5306524.215) ( 44d13'51.88"E, 42d57'51.12"N)
Band 1 Block=1656x4 Type=Byte, ColorInterp=Gray

So the question is:

How to re-project wide-angle camera? Looks like gdal incorrectly process it. I.e. pixel size in the middle of the picture is ~1km. While gdalinfo makes it ~2km. Which might be ok on the edges.

Another question:

Why only “-tps” can handle non-linear distortions on the edges?

One Answer

Issue was on my side. When I was generating GCPs, I took (originalImageHeight / someRate) points. Effectively missing the last line. So the correct number of lines should be (originalImageHeight / someRate) + 1.

enter image description here

Answered by Andrey Rodionov on December 4, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP