How to show long and lat values on a map?
I have a projection file EPSG:3410
that I want to simply plot using R. Since the projection is in meters I have a problem with wide and long values (axes). You can see on the map that instead of 90 -90 180 -180 there are other numbers.
the code to build this map:
con1 <- file("C:\\Users\\data.bin","rb")
r = raster(y)
extent(r) = extent(c(xmn=-17334194,xmx=17334194,ymn=-7356860,ymx=7310585))
plot(r)
plot(wbuf, add=TRUE)
- How can I replace these unusual values in the map from 90 -90 180 -180 (so no changes to the original file)
- How can I get rid of the space at the top and bottom of the map?
source to share
Here are some workarounds to avoid whitespace and add axis labels. Since both cards are in meters, it is not easy to construct lat long marks.
library(raster) # for xmin, ymin etc.
# plot raster without axes and box
plot(r, asp=1, axes=F, box=F)
# add polygons
plot(wbuf, add=TRUE, axes=F)
# draw bounding box at raster limits
rect(xmin(r), ymin(r), xmax(r), ymax(r))
To add lat-lon ticks, we need to convert from EPSG: 3410 (meters) projection to lat-lon (degrees). We can create dummy spatial points in the latlon and transform them to 3410 and then extract the transformed coordinates to add ticks
# y axis:
x = -180
y = seq(-90,90,30) # you can adapt this sequence if you want more or less labels
S <- SpatialPoints(cbind(x,y), proj4string = CRS("+proj=longlat +datum=WGS84"))
S2<- spTransform(S, CRS(" +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +no_defs "))
axis(side = 2, pos=xmin(r), at = S2@coords[,'y'], lab=y, las=1)
# x axis
y = -90
x = seq(-180,180,30)
S <- SpatialPoints(cbind(x,y), proj4string = CRS("+proj=longlat +datum=WGS84"))
S2<- spTransform(S, CRS(" +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +no_defs "))
axis(side = 1, pos=ymin(r), at = S2@coords[,'x'], lab=x)
Another option is to reprogram the raster so that the maps are lat-lon (degrees):
# reproject and plot raster
r2 = projectRaster(r, crs=CRS("+proj=longlat +datum=WGS84"))
plot(r2, axes=F, box=F)
# there an option to disable internal boundaries, which is nice
map(database = 'world', regions = '', interior = F, add=T)
# still have to add bounding box and axes
rect(xmin(r), ymin(r), xmax(r), ymax(r))
axis(1, pos = ymin(r2))
axis(2, pos= xmin(r2), las=1)
source to share