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.rowling...@gmail.com [mailto:b.rowling...@gmail.com] De la part de Barry Rowlingson Envoyé : mardi 1 octobre 2013 17:43 À : Amelie LESCROEL Cc : r-sig-geo@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.lescr...@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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
