[R-sig-Geo] Import a GeoTIFF file and plot it on a map with lat-long coordinates

Amelie LESCROEL Amelie.LESCROEL at cefe.cnrs.fr
Tue Oct 1 15:26:05 CEST 2013


Dear all,
I'm having trouble creating a map with lat-long coordinates from a (supposedly?) GeoTIFF file and I would gratefully welcome any suggestion of more advanced users. This file is a satellite image of the southern Ross Sea, Antarctica, which can be found here: http://lance-modis.eosdis.nasa.gov/imagery/subsets/?project=antarctica&subset=RossSea.2004356.aqua.500m
As proposed on this page, I downloaded the .tif file as well as the .prj file, then tried to import and plot it in R.

> library(raster)
Le chargement a nécessité le package : sp
> library(rgdal)
rgdal: version: 0.8-11, (SVN revision 479M)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files: C:/Users/lescroel/Documents/R/win-library/2.15/rgdal/gdal
GDAL does not use iconv for recoding strings.
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files: C:/Users/lescroel/Documents/R/win-library/2.15/rgdal/proj
> setwd(./.)
> sat <- raster("RossSea_2004356_aqua_500m.tif")
> sat
class       : RasterLayer 
band        : 1  (of  3  bands)
dimensions  : 1800, 1400, 2520000  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 1400, 0, 1800  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : F:\Mes Documents\Postdoc\Articles\Phenotypic plasticity of foraging efficiency\PLOS ONE\RossSea_2004356_aqua_500m.tif 
names       : RossSea_2004356_aqua_500m 
values      : 0, 255  (min, max)
> plot(sat, col=grey(0:255/255))

It works but all the projection info is lost. And when I plot it, I end up with axes going from 0 to 1400 for x and 0 to 1800 for y (as specified in the extent of the sat object), not lat-long coordinates.
I tried reading the .tif file with readGDAL but I get an error message:

> sat2<-readGDAL("C:/Users/lescroel/Desktop/RossSea.2004356.aqua.500m.tif")
C:/Users/lescroel/Desktop/RossSea.2004356.aqua.500m.tif has GDAL driver GTiff 
and has 1800 rows and 1400 columns
Message d'avis :
In readGDAL("C:/Users/lescroel/Desktop/RossSea.2004356.aqua.500m.tif") :
  GeoTransform values not available

I tried specifying a projection from what I could see of the .prj file (PROJCS["Lambert Azimuthal Sphere",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-170.0],PARAMETER["Latitude_Of_Origin",-75.0],UNIT["Meter",1.0]]), but I'm not even sure that I "translated" it correctly:

> sat.geo <- CRS("+proj=laea +lat_0=-75 +lon_0=-170 +a=6371000 +units=m +no_defs")
> proj4string(sat) <- sat.geo
> sat
class       : RasterLayer 
band        : 1  (of  3  bands)
dimensions  : 1800, 1400, 2520000  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 1400, 0, 1800  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-75 +lon_0=-170 +a=6371000 +units=m +no_defs 
data source : F:\Mes Documents\Postdoc\Articles\Phenotypic plasticity of foraging efficiency\PLOS ONE\RossSea_2004356_aqua_500m.tif 
names       : RossSea_2004356_aqua_500m 
values      : 0, 255  (min, max)

Is there a way to read the .prj directly from R? How should I do to be able to map this image with long-lat coordinates? I guess that lots of the required information to do this could be gathered from the metadata of this file (http://lance-modis.eosdis.nasa.gov/imagery/subsets/?project=antarctica&subset=RossSea.2004356.aqua.500m.met) but I'm not sure how to implement this and all my previous attempts failed.

Thanks in advance for your help and do not hesitate to direct me towards an online resource / reading material that I could have overlooked.

Best regards,

Amélie

> sessionInfo()
R version 2.15.3 (2013-03-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252    LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_0.8-11  raster_2.1-49 sp_1.0-13    

loaded via a namespace (and not attached):
[1] grid_2.15.3     lattice_0.20-13 tools_2.15.3
_______________________________________________ 
Amélie Lescroël 
Centre d'Ecologie Fonctionnelle et Evolutive 
CNRS - UMR 5175 
1919 Route de Mende 
F-34293 Montpellier Cedex 05 



More information about the R-sig-Geo mailing list