Tuesday 21 January 2020

raster - How can I fix a geoTiff with an incorrect CRS?


Trying to load and plot a raster (it's a DEM) in R, but having problems.


Using R I read in a raster file:


library(raster)

dem <- raster("DEM_Victoria.tif") # >70 MB

print(dem)

class : RasterLayer
band : 1 (of 3 bands)
dimensions : 5444, 6007, 32702108 (nrow, ncol, ncell)
resolution : 100, 100 (x, y)
extent : 176615.6, 777315.6, 5650397, 6194797 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +no_defs
data source : DEM_Victoria.tif
names : DEM_Victoria

values : 0, 255 (min, max)

My confusion here is the extent vs the CRS. The extent shows coordinates from a UTM Zone 55S system, while the CRS reports it as longlat.


Similarly, gdalinfo reports:


Driver: GTiff/GeoTIFF
Files: DEM_Victoria.tif
Size is 6007, 5444
Coordinate System is:
GEOGCS["GDA94",
DATUM["unknown",

SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (176615.587999997427687,6194797.138000000268221)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND

Corner Coordinates:
Upper Left ( 176615.588, 6194797.138) (Invalid angle,Invalid angle)
Lower Left ( 176615.588, 5650397.138) (Invalid angle,Invalid angle)
Upper Right ( 777315.588, 6194797.138) (Invalid angle,Invalid angle)
Lower Right ( 777315.588, 5650397.138) (Invalid angle,Invalid angle)
Center ( 476965.588, 5922597.138) (Invalid angle,Invalid angle)
Band 1 Block=6007x1 Type=Byte, ColorInterp=Red
Band 2 Block=6007x1 Type=Byte, ColorInterp=Green
Band 3 Block=6007x1 Type=Byte, ColorInterp=Blue


Is this real? A problem? (well it is to me, but perhaps because I don't understand something) Fixable?


I tried 'Assign projection' in QGIS, which uses this command:


gdalwarp -t_srs EPSG:28355 -srcnodata -32768  DEM_Victoria.tif DEM_Victoria.tif.tmp

but it failed with:



ERROR 1: Too many points (12100 out of 12100) failed to transform, unable to compute output bounds.




Answer



If I understand this correctly, you're not translating the coordinates, you need to assign the correct CRS. gdal_translate can do this:



gdal_translate -a_srs EPSG:28355 DEM_Victoria.tif DEM_Victoria.tif.tmp

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...