7

For several months, I have been georeferencing files manually using the Georeferencer plugin within QGIS. I am now attempting to automate this process using GDAL. The Georeferencer provides a "Generate GDAL script" option, which creates the following commands:

gdal_translate -of GTiff -gcp 774.73 -338 -88.52 34.37 -gcp 1012.17 -215.136 -81.22 27.48 -gcp 972.709 -486.055 -77.76 36.67 asciiMAE20121001.asc testoutput1
gdalwarp -r near -order 1 -co COMPRESS=NONE testoutput1 testoutput2

When I run this script, however, the georeferenced raster is placed about 40 degrees further west than where it should be, and 7 degrees further south. Manually using the Georeferencer plugin continues to georeference the raster correctly. I think the problem is with the gdalwarp, because the raster is just shifted too far from where it should be. The placement of the values within it appears to be correct, so it looks like the gdal_translate worked correctly.

Does anyone know why the GDAL script produced by the Georeferencer would produce a raster displaced from the raster produced manually by the plugin?

AndreJ
  • 76,698
  • 5
  • 86
  • 162
Chris S.
  • 101
  • 1
  • 4
  • Take gdalinfo from both variants. It may be that the command line commands do not set the projection correctly. – user30184 Oct 10 '14 at 14:53
  • 1
    btw, I recently asked a question using the same process. Might be a useful comparison: http://gis.stackexchange.com/questions/116672/georeferencing-a-raster-using-gdal-and-python – djq Oct 10 '14 at 14:59
  • Also, another useful example: https://www.mapbox.com/blog/georeferencing-satellite-geoeye/ – djq Oct 10 '14 at 14:59
  • 1
    Have a try by adding -s_srs epsg:1234 -t_srs epsg:1234into your generated gdalwarp command. Change 1234 into the projection of your ground control points. – user30184 Oct 10 '14 at 15:01
  • I added the above -s_srs and -t_srs to epsg:4326, which should match the WGS 84 that I use with the manual georeference. Taking the gdalinfo of both files, I see that the manually created file has a coordinate system of "'", rather than the WGS 84 specified. What does the "'" specify? – Chris S. Oct 10 '14 at 15:33
  • The character before the single quote in my above comment should be a grave accent; it is not a space like displayed above. – Chris S. Oct 10 '14 at 15:45
  • If you use -s_srs and -t_srs you should get image that is tagged to be in EPSG:4326 or otherwise something has failed. This is how the beginning of projection info looks for me Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], – user30184 Oct 10 '14 at 16:01
  • Yes, the script appears to be tagging it correctly, but the manual georeference has the strange coordinate system. When doing the manual georeference, I specify WGS 84 in the QGIS GUI, and I see it when I check the raster's properties. I'm not sure why the gdalinfo statement does not show this coordinate system when I look at the manually-created raster. – Chris S. Oct 10 '14 at 16:32
  • Obviously gdalinfo does not show the coordinate system because QGIS does not save that info into geotiff. It may be that QGIS, when it opens the warped file automatically into the project, knows where to place it even if the projection is unknown because it is opened from the georeferencing task. Does the file which you save properly into EPSG:4326 open in a correct place? – user30184 Oct 10 '14 at 16:42
  • No, the file created by the warp command opens much further west (and slightly further south), even now that I have specified EPSG:4326. The manually created raster, as well as all of my shapefiles, are 40 degrees longitude to the east of the raster I create through the warp command (and slightly further north too). – Chris S. Oct 10 '14 at 16:53
  • Can you share the .asc source file via dropbox or similar? – AndreJ Oct 11 '14 at 04:38
  • 2
    Yes, we would need test data for repeating the issue. What I think is that QGIS writes wrong pixel/line offsets to gdalwarp command. Notice that in all your GCPs the line offsets are negative which means for GDAL that they are north to the top left pixel, thus totally outside the image area. QGIS is using negative offset natively in the top-down direction but it usually turns the sign correctly for the GDAL commands. However, I noticed that if original image is already georeferenced the QGIS behaves all wrong - it takes georeferenced coordinates from the source and not pixels/lines. – user30184 Oct 11 '14 at 10:41
  • Thank you, I believe your analysis is correct, which means I probably can't remedy this problem through my warp command. With this in mind, I found another way to tackle the problem. I used the python script from How to translate (reposition) a .tif raster layer? to shift the raster to the correct lat/long, and now it matches the manually created raster. – Chris S. Oct 11 '14 at 14:50
  • The strange thing for me is when I georeference a jpg in QGIS, the srcY in the GCP table are all negative but in the GDAL commandline they are inverted to positive, giving the right referencing in the end. Seems like a bug if an .asc file is loaded. – AndreJ Oct 13 '14 at 19:05
  • You have negative pixel lines. Shoud't line and pixel be positive, alays? – wingnut May 13 '21 at 09:25

3 Answers3

2

Assuming you have saved the output of your GCP points, can you compare the values in this file with the auto-generated script? I've noticed that the script rounds values; I wonder if this could be introducing a source of error.

djq
  • 16,297
  • 31
  • 110
  • 182
  • The script does indeed round the values, but it rounds to the nearest hundredth, so the rounded values are very close of those of the saved file. Just in case, I swapped in the original points, but there was not a change. Thank you for the suggestion though! – Chris S. Oct 10 '14 at 15:27
1

I have encountered this problem recently as well so I checked the pixel coordinate system used by the QGIS Georeferencer and noticed that the origin (upper left corner) is at (1000, 1000) instead of (0, 0).

You might want to check if your origin is not (0, 0) and apply necessary offsets to the pixel coordinates and see if gives you the right output (it worked for me).

R.K.
  • 17,405
  • 3
  • 59
  • 110
0

The comment by @user30184 indicates the problem. The manual shows that the arguments to gcp are "pixel line easting northing elevation". For the qgis generated gdal script to be correct, you must have loaded a raster without a reference system. Otherwise the first two arguments for each -gcp will be wrong.

I removed the spatial information from my tif with this:

gdal_translate -co PROFILE=BASELINE in.tif out.tif

Then I loaded out.tif and found the ground control points using qgis and exported the script.

Alternatively you could convert the first two arguments of each of your existing gcps to the pixel and line of the image.

Tedward
  • 101
  • 1