2

I have images georeferenced to EPSG:32612 and coordinates in EPSG:4236. I see that both use the WGS84 ellipsoid. I want to convert the images to EPSG:4326 to mark positions in the images using lat/long, so I used gdalwarp to warp the image, like so:

gdalwarp -s_srs EPSG:32612 -t_srs EPSG:4326 infile.tif outfile.tif

I expected this to simply change the metadata in the TIFF because, so far as I understand, these two coordinate systems are only a translation and scale apart.

I then imported both infile.tif and outfile.tif into QGIS to double-check, and was surprised to see that the warp had actually changed the pixel values. Below you can see me swapping between the two images with a Google Satellite layer in the background.

comparing in QGIS

So, is my understanding incorrect? Is a reprojection using the same ellipsoid not an identity operation on pixel values? If not, why not? If it is, why does QGIS give noticably different pixel values when rendering each image? Is gdalwarp simply not pixel-perfect? Or is QGIS doing a bad "nearest" reprojection for display?

Practically speaking, if it should be the identity, is there a way to modify just the meta-data for a tif image using GDAL, or rasterio? (I tried rasterio.warp.calculate_default_transform, but it requires RPCs and anyway ends up more than 10 metres incorrect)

Note: I know that I could translate the coordinates into EPSG:32612. But, this effect would still occur and I'm just trying to understand why the positions are slightly different.

3
  • 3
    These two coordinate systems are not only a translation and scale apart. One is Cartesian 2D CS. Axes: easting, northing (E,N). Orientations: east, north. UoM: m, the other Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree. In reprojection the pixels must be resampled.
    – user30184
    Aug 16, 2021 at 11:43
  • 3
    UTM is a projection from the GCS. While they share an ellipsoid, GCS WGS84 uses angular units (degrees) and UTM Zone 12N uses meters. Converting between UTM 12N and UTM 12S would be a case of translation.
    – Vince
    Aug 16, 2021 at 12:07
  • 1
    No need for reprojectio if you only want to get lat/lon values. Load the raster in EPSG:32612, than right-click on the canvas to get coordinates. In newer QGIS versions, you can choose in what CRS you get the coordinate values. If you load geometries in 4326, QGIS will place them correctly. Even if you did not explain your workflow, I strongly suppose you don't nee any re-projection.
    – Babel
    Aug 16, 2021 at 12:07

2 Answers 2

4

Your basic error is that you think a rectangular UTM image would easily reproject to an identical WGS84 rectangle. But latitudes and longitudes are bended in any UTM projection, and vice versa a kilometer grid in a WGS84 projection. Same goes for the extent of every single pixel.

So GDAL has to create a new rectangle that is slightly expanded to get all your pixels into. And it has to resample the content of the UTM file pixels to get values for the WGS grid.

If your source file has paletted colours, this goes totally berserk. You might get better values if you transform to a RGB image before reprojecting.

Regarding points, your understanding of identity is however correct. For lines and polygons, you will see differences if you densify the geometry before or after the reprojection.

For illustration, I have created a 10m-grid in UTM12 (blue) over Salt Lake City, and put a WGS84 grid in red with the same cell size on top of it: enter image description here

You see that they align only in the left bottom of the map. Each square should represent a pixel of a low-res satellite image.

6
  • Thanks for your answer. Just checking my understanding. Every pixel value in a satellite image is the average colour of a small area on the surface of the Earth. When you reproject from EPSG:32612 (metres) to EPSG:4326 (degrees) the area that each pixel in the image represents changes shape slightly, which means that a resampling is required. The real world point described by the origin of each pixel in the image is unchanged between these projections, but the box of area in the real world that the pixel data represents has changed. Is that correct? Aug 17, 2021 at 5:25
  • Almost. Gdal spans a completely new grid of pixels over your original image. The origin (better: center) of the pixels don't have to be on the same spot. The pixel value is then calculated from the values of the underlying original pixels. You can define the new pixel size yourself with parameters. So you can create a one-degree-grid interpolating from a one-meter UTM grid. This will obviously look different from the original image. GDAL tries to preserve the number of pixels in height and length as much as possible.
    – AndreJ
    Aug 17, 2021 at 6:11
  • Actually, I am fairly lost. Why must it lay a completely new grid? Imagine projecting a grid of pixels from the satellite's image plane onto the surface of a perfect ellipsoid. How we choose to measure the locations of the pixel centers on the surface doesn't affect the projection at all, right? That is, the areas which each pixel collected light from are still the same physical areas regardless of whether we measure the locations in degrees or meters. ... Wait, are these images assumed to be orthorectified? Aug 17, 2021 at 6:57
  • I updated my answer.
    – AndreJ
    Aug 17, 2021 at 12:34
  • Thanks for your effort here; you have clarified a lot. I apologise, but I am still not sure I fully understand the why. I understand that there is no uniform grid that is the same in both CRS. And if I have an image where each pixel covers a uniform area in one CRS, and I want to preserve this property in the new CRS, I can understand the need to resample. But if my image hasn't been aligned to either uniform grid yet (i.e. hasn't been orthorectified), I don't know why it should resample. I feel like it shouldn't. Aug 18, 2021 at 2:37
0

The resampling is to maintain orthorectification between the two CRS.

If the images are not orthorectified, then measuring the image boundaries in different units doesn't affect the pixel information. However, reprojection typically implies using orthorectified images and going from an image orthorectified on EPSG:4326 to an image orthorectified on EPSG:32612 involves a resample operation.

Why?

EPSG:4326 is measured in degrees, but EPSG:32612 is measured in metres. EPSG:4326 are polar coordinates, but EPSG:32612 are cartesian coordinates. gdalwarp assumes that you have provided an orthorectified image. That is, each pixel represents some unit area. e.g. each pixel is 0.001 degrees across for EPSG:4326 vs each pixel is 1m across in EPSG:32612.

Imagine laying a uniform grid on the surface of the Earth for each of these CRS. Clearly these grids will never overlap perfectly, regardless of your pixel size (because of cartesian vs polar). For gdalwarp to produce an orthorectified image in the new CRS, it must resample the image on a new grid. It will choose the pixel size such that there is minimal difference and the resampled image looks most similar. But the area covered by each pixel has changed.


The ellipsoid only describes the surface that these projections are coming from. If the ellipsoid changed, that would be an extra reason to resample the image.

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.