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

[Amelie LESCROEL]

Thanks a lot for your solution, which works perfectly well, Barry. I should have worked with the .jpg rather than the .tif file to be able to read the projection and explore the function projectRaster(). Now, I just have to fiddle with the graphical parameters so that the y axis is not too far on the left from the image border.
Amélie

-----Message d'origine-----
De : b.rowlingson at gmail.com [mailto:b.rowlingson at gmail.com] De la part de Barry Rowlingson
Envoyé : mardi 1 octobre 2013 17:43
À : Amelie LESCROEL
Cc : r-sig-geo at r-project.org
Objet : Re: [R-sig-Geo] Import a GeoTIFF file and plot it on a map with lat-long coordinates

On Tue, Oct 1, 2013 at 2:26 PM, Amelie LESCROEL <Amelie.LESCROEL at cefe.cnrs.fr> wrote:


> 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.

 If you get the "Download JPG image with ancillary files (.zip)" and unzip it as well as the tif you'll get a .jpg, a .jgw, and a .prj file. Gdal will use those to give the JPG file a coordinate system:

 > j=stack("./RossSea.2004356.aqua.500m.jpg")
 > projection(j)
[1] "+proj=laea +lat_0=-75 +lon_0=-170 +x_0=0 +y_0=0 +a=6371000
+b=6371000 +units=m +no_defs"

 but as you noticed this doesn't work for GeoTIFFs:

> t=stack("./RossSea.2004356.aqua.500m.tif")
> projection(t)
[1] "NA"

never mind, as long as they have the same projection, we can just:

projection(t)=projection(j)
plotRGB(t)

and we don't need 'j' any more.

There's probably a proper way to do this, but this works.

Note that converting to long-lat is going to need a warp operation. I normally do this in two steps:

1. compute the extent in the new coords:

e=projectExtent(t,CRS("+init=epsg:4326"))

2. do it:

tl=projectRaster(t,e,method="ngb")

[ngb is a bit quicker than the default, read the docs and decide]

plotRGB(tl)

should then have it in lat-long. Its not much of a stretch since its a smallish area.

Barry