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