[R-sig-Geo] drawing a vector map on top of an aerial photograph

[Roger Bivand]

On Tue, 1 Apr 2008, Bjarke Christensen wrote:

> Hi,
>
> I have a bunch of 256x256 pixel aerial photographs that span a territory
> for which I also have simple arcview vector maps. I'd like to be able to
> plot these vector maps along with some vector shapes that result from my
> calculations, and to use the aerial photos as background.
>
> Roger Bivand has answered a similar question on
> http://tolstoy.newcastle.edu.au/R/e2/help/06/11/6178.html - but I can't
> seem to get the solution to work. Specifically,
> image(readGDAL("image.jpg")) produces an orange-tinted image which
> naturally has no coordinates (since it's just a jpeg).

It is plotting the first band using default heat.colors(), hence the 
orange appearance. If you want three bands, you'll need to add red=, 
green= and blue= arguments to the call to image().

>
> So my problem is firstly, how to attach coordinates to the raster image
> when loading it,

Can you use a format that has the correct coordinates in it already, such 
as Geotiff? Do you know the correct coordinates for the jpeg images? Do 
you know how they are represented in terms of the coordinate reference 
system? If so, something might be possible. Do you, for example have a 
World file for the jpeg, and is it in the same directory - if it was, GDAL 
would use it, see http://www.gdal.org/frmt_jpeg.html. See:

http://www.nps.gov/gis/av3_online/documents/section7/docs/section7_pg3.html

for how you might create one.

> secondly, how to resample the image up to the same scale
> as the vector map,

Doing image() of the correctly configured raster data (with both known 
coordinates and coordinate reference system) makes this easy - just do 
plot(..., add=TRUE) to plot the vector data (providing the coordinate 
reference systems are the same, otherwise use spTransfrom() in rgdal).

> and thirdly, how to 'glue' the tiles together to span
> the area being mapped.

If this is just for display, why not just set xlim= and ylim= in the first 
image(), then do image(..., add=TRUE) to fill in. There are rbind() 
methods for SpatialPixelsDataFrame objects in sp, so you could cast to 
those, or straight to SpatialPointsDataFrame objects, which can then be 
cast back to a SpatialGridDataFrame object when done.

Roger

>
> Any help would be greatly appreciated!
> Bjarke.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no