3

I have a satellite image as a PNG with known coordinates (bottom, left, top, right). I want to draw a shapefile on it (the limits of a county) in R.

My first thought was to convert the PNG (it's read as a three-layered matrix, for red, green and blue) to a raster. But creating a three-layered raster may be an overkill.

And simpler solution to do this?

I found the function rasterImage(), and I'm doing like that:

Map0 <- readPNG('myMap0.png')
cM <- c(xmin,ymin,xmax,ymax)
png('myMap.png',10240,8610)
par(mar=c(0,0,0,0))
plot(myShapeFile,border=NA,xlim=c(cM[1],cM[3]),ylim=c(cM[2],cM[4]))
rasterImage(Map0,cM[1],cM[2],cM[3],cM[4])
plot(myShapeFile,add=T,border='white',lwd=5)
dev.off()

It's almost fine, except for a white border around the image. Since I used the dimensions of the original image (10240x8610), this border is probably distorting the new image, because the painted area is smaller than the original. I had this problem before, but could never solve it.

I thought the problem was with png(), so I did the following:

par(mar=c(0,0,0,0),bg=NA)
plot(myShapeFile,border='white',xlim=c(cM[1],cM[3]),ylim=c(cM[2],cM[4]),lwd=5)
dev.print(width=10240,height=8610,'myBorder.png',dev=png)

So I could try to merge the two images as layers in Gimp, but I see that the border remains, causing the shapefile to be drawn in the wrong position (it aligns perfectly with myMap.png, but not with myMap0.png). In Google Earth I know that the coordinates in cM are precise, so, since I'm giving mar=0, and xlim and ylim should refer to the limits of the plot, why this margin keeps on disturbing the image created?

Following @mdsumner comment, I tried the following:

map0 <- brick('myMap0.png')
extent(map0) <- c(cM[1],cM[3],cM[2],cM[4])
png('myMap.png',10240,8610)
par(mar=c(0,0,0,0))
plotRGB(map0,maxpixels=10240*8610)
plot(myShapeFile,add=T,border='white',lwd=5)
dev.off()

Now it's almost right, except for the left and right white margins, each 27 pixels wide).

5
  • 1
    @Aaron Please see my edit.
    – Rodrigo
    Aug 30, 2017 at 1:58
  • Try raster package. Map0 <- brick('myMap0.png'); extent(Map0) <- c(xmin, xmax, ymin, ymax); plotRGB(Map0); plot(myShapeFile, add = TRUE)
    – mdsumner
    Aug 30, 2017 at 3:17
  • @mdsumner Still no success. Please see my edit.
    – Rodrigo
    Aug 30, 2017 at 3:34
  • I reckon your extents are wrong, or your shapefile is in a different CRS. Fudging map data in a PNG with free floating numbers to georeference it from scratch is a very tenuous situation, and to work through possible issues will require a lot more information. At minimum, include the plot/s you see, and print summaries of the shapefile object and the raster object - the ones you plot with.
    – mdsumner
    Aug 30, 2017 at 15:14
  • @mdsumner You' right. I have re-made the calculations, and there was a small error in the bottom coordinate (since it's a mosaic of Google Maps imagery, I had to remove the 50px of the logo, and forgot to update the coordinate accordingly). You may want to include your suggestion as the accepted answer. Thank you!
    – Rodrigo
    Aug 30, 2017 at 15:48

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.