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

[Amelie LESCROEL]

Actually, now that I'm trying (unsuccessfully) to add a SpatialPolygon to my plotted raster, I realize that something went wrong in the reading or projection process of the .tif file, which is why you said " since its a smallish area". Well, it's not: it should go from about 70 to 80°S and from 160°E to 170°W. But when I do plotRGB(tl, axes=TRUE), the axes go only from -75.00000 to -74.98381 and from -170.0000 to -169.9514. Yikes!
I guess there is something wrong with the extent: 

> extent(t)
class       : Extent 
xmin        : 0 
xmax        : 1400 
ymin        : 0 
ymax        : 1800 
> extent(j)
class       : Extent 
xmin        : -948250 
xmax        : -248250 
ymin        : -618000 
ymax        : 282000 
> extent(tl)
class       : Extent 
xmin        : -170 
xmax        : -169.9514 
ymin        : -75 
ymax        : -74.98381

The metadata of the .tif file say:

projection center lon: -170.0000
projection center lat: -75.0000
image center lon: +168
image center lat: -75.5
UL lon: +163.5539
UL lat: -70.5374
UR lon: -177.3704
UR lat: -72.3252
LR lon: +176.6394
LR lat: -80.3056
LL lon: +147.4332
LL lat: -77.3404
UL easting (km): -948.2500
UL northing (km): +282.0000

Any hint about how I could get this right?

Thanks for your consideration,

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