3

I'm trying to reproject a EPSG 3035 geotif to a EPSG 3857 using gdalwarp as follows:

gdalwarp -tap -tr 100 100 -s_srs EPSG:3035 -t_srs EPSG:3857 -co 
SPARSE_OK=TRUE -co predictor=2 -of GTiff -co COMPRESS=DEFLATE -co 
TILED=YES -multi -wo NUM_THREADS=ALL_CPUS source.tif target.tif

The source tif is generated in GRASS (7.2) and correctly displayed in qgis while the target one contains only nan. The source tif has the following metadata (gdalinfo):

PROJCS["ETRS89 / LAEA Europe",
GEOGCS["grs80",
    DATUM["European_Terrestrial_Reference_System_1989",
        SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6258"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",52],
PARAMETER["longitude_of_center",10],
PARAMETER["false_easting",4321000],
PARAMETER["false_northing",3210000],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (944300.000000000000000,5413400.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.2.1 with GDAL 2.1.2
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  944300.000, 5413400.000) ( 54d53' 7.79"W, 56d24'16.47"N)
Lower Left  (  944300.000,  942200.000) ( 23d33' 5.25"W, 24d48'30.69"N)
Upper Right ( 6525800.000, 5413400.000) ( 60d11'18.71"E, 64d12'38.48"N)
Lower Right ( 6525800.000,  942200.000) ( 32d24'30.12"E, 28d34'52.61"N)
Center      ( 3735050.000, 3177800.000) (  1d33'27.09"E, 51d24'20.50"N)
Band 1 Block=55815x1 Type=Float32, ColorInterp=Gray
  NoData Value=nan

In order to generate the source tif I used some other tifs that show the following metadata:

PROJCS["ETRS89 / LAEA Europe",
GEOGCS["ETRS89",
    DATUM["European_Terrestrial_Reference_System_1989",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6258"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4258"]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",52],
PARAMETER["longitude_of_center",10],
PARAMETER["false_easting",4321000],
PARAMETER["false_northing",3210000],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3035"]]
Origin = (944300.000000000000000,5413400.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  944300.000, 5413400.000) ( 54d53' 7.79"W, 56d24'16.47"N)
Lower Left  (  944300.000,  942200.000) ( 23d33' 5.25"W, 24d48'30.69"N)
Upper Right ( 6525800.000, 5413400.000) ( 60d11'18.71"E, 64d12'38.48"N)
Lower Right ( 6525800.000,  942200.000) ( 32d24'30.12"E, 28d34'52.61"N)
Center      ( 3735050.000, 3177800.000) (  1d33'27.09"E, 51d24'20.50"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=0

They are similar but not equals... also if the proj string is the same:

+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Why can I not properly convert the source tif?

Are there some bugs in GRASS or gdal?

For complete information this is the output of gdainfo on the target tif:

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 = (-6109900.000000000000000,11708100.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.2.1 with GDAL 2.1.2
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-6109900.000,11708100.000) ( 54d53'10.20"W, 71d52'27.93"N)
Lower Left  (-6109900.000, 2852200.000) ( 54d53'10.20"W, 24d48'29.38"N)
Upper Right ( 6700300.000,11708100.000) ( 60d11'23.35"E, 71d52'27.93"N)
Lower Right ( 6700300.000, 2852200.000) ( 60d11'23.35"E, 24d48'29.38"N)
Center      (  295200.000, 7280150.000) (  2d39' 6.58"E, 54d34'36.27"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=nan

UPDATE: using --debug on and suppressing -s_srs. The resulting tif is the same:

GDAL: GDALOpen(results_last_update/prediction_14.tif, 
this=0x750de5a6d0) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGRCT: Source: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 
+y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Target: +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
OGRCT: Source: +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
OGRCT: Target: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 
+y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000         +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Target: +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
OGRCT: Source: +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
OGRCT: Target: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
Creating output file that is 128102P x 88559L.
GDAL:     GDALDriver::Create(GTiff,idi/prediction_14_2.tif,128102,88559,1,Float32,0x750de30b70)
Processing input file results_last_update/prediction_14.tif.
WARP: Copying metadata from first source to destination dataset
GDAL: GDALDefaultOverviews::OverviewScan()
Using internal nodata values (e.g. nan) for image     results_last_update/prediction_14.tif.
Copying nodata values from source     results_last_update/prediction_14.tif to destination     idi/prediction_14_2.tif.
WARP: srcNoData=nan dstNoData=nan
WARP: calling GDALSetRasterNoDataValue() for band#0
OGRCT: Source: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Target: +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
OGRCT: Source: +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
OGRCT: Target: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Target: +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
OGRCT: Target: +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
OGRCT: Source: +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
OGRCT: Target: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +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
OGRCT: Target: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
GDAL: Start chunk 0.
GDAL: Start chunk 1.
0GDAL: GDAL_CACHEMAX = 797 MB
GDAL: Finished chunk 0.
GDAL: Start chunk 2.
...

I have tried with a smaller area and gdalwarp seems to works. Maybe the boundin box of the resulting tif is not correct? I'm trying to find the largest bounding box.

  • If I use gdalwarp on an input tif, this is correctly reprojected and displayed so I think the problem is caused by GRASS (r.gdal.out ) – Andrea Gobbi Sep 20 '17 at 12:41
  • Have you tried to use gdalwarp without -s_srs ? Adding --debug on at the end of the command may give more info. – user30184 Sep 20 '17 at 12:51
  • have you tried proj string instead of epsg ? – radouxju Sep 20 '17 at 14:07
  • gdalwarp -tr 100 100 -tap -t_srs "+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" -co SPARSE_OK=TRUE -co predictor=2 -of GTiff -co COMPRESS=LZW -co TILED=YES -multi -wo NUM_THREADS=ALL_CPUS source.tif target.tif same result – Andrea Gobbi Sep 20 '17 at 14:41
  • Your original tifs had a nodata value of zero. Perhaps you should assign that when creating your source tif. – AndreJ Sep 20 '17 at 17:55
1

Nevermind. At the end I user R to read the target tif. I noticed that there was values. It is just a QGIS problem of visualization. Thanks anyway

  • For such type of pixel value you need to create statistics (min, max, std deviations etc.) using QGIS or gdalinfo -stats command. And use min-max histogram stretch. By default QGIS not stretch you Float values to visible 0 - 255 RGB values. Anyhow you need to setup raster visualisation style. – Dmitry Baryshnikov Sep 21 '17 at 9:21

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.