Gdal_calc outputs image with no projection

Question

I am trying to convert a pixel value above 0.02 to be nodata. However, the output file comes out with no projection. When I assign projection, QGIS loads it in the wrong place.

gdal_calc.py -A "000-IMG_0037_4_gcp.tif" --outfile="000-IMG_0037_4_gcp_noglint8.tif" --calc="((A<=0.02)*(A>0))*A+((A>0.02)*0)" --NoDataValue=0

The gdalinfo of the original file:

Files: /drive/000-IMG_0037_4_gcp.tif
Size is 1048, 1422
Coordinate System is `'
GCP Projection = 
PROJCS["WGS 84 / UTM zone 16N",
    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["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32616"]]
GCP[  0]: Id=1, Info=
      (0,0) -> (280191.450981,3352631.04297,0)
GCP[  1]: Id=2, Info=
      (1048.9984,0) -> (280298.687984,3352631.04297,0)
GCP[  2]: Id=3, Info=
      (1048.9984,1422.4381) -> (280298.687984,3352484.8363,0)
GCP[  3]: Id=4, Info=
      (0,0) -> (280191.450981,3352631.04297,0)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1422.0)
Upper Right ( 1048.0,    0.0)
Lower Right ( 1048.0, 1422.0)
Center      (  524.0,  711.0)
Band 1 Block=1048x1 Type=Float32, ColorInterp=Gray

The output file gdalinfo:

Driver: GTiff/GeoTIFF
Files: /drive/000-IMG_0037_4_gcp_noglint8.tif
Size is 1048, 1422
Coordinate System is `'
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1422.0)
Upper Right ( 1048.0,    0.0)
Lower Right ( 1048.0, 1422.0)
Center      (  524.0,  711.0)
Band 1 Block=1048x1 Type=Float32, ColorInterp=Gray
  NoData Value=0

Is there a way to make sure that the output image comes with the same srs of the input raster?

Answer 1

Your original file isn't actually georeferenced, it uses GCPs. I'm not sure why gdal_calc ignores the GCPs, but you can workaround using gdalwarp to georeference the original file.

The following should work:

gdalwarp 000-IMG_0037_4_gcp.tif 000-IMG_0037_4_warped.tif
gdal_calc.py -A "000-IMG_0037_4_warped.tif" --outfile="000-IMG_0037_4_warped_noglint8.tif" --calc="((A<=0.02)*(A>0))*A+((A>0.02)*0)" --NoDataValue=0