2

I have had trouble mapping categorical points on to a map using rworldmaps or the maps package in R.

The original data was Aragonite saturation levels in the ocean on an .adf file which was in Lambert Cylindrical Equal Area format ("+proj=cea +datum=wgs84 +lon_0=-160.0 +lat_ts=0.0 +x_0=0.0 +y_0=0.0"). When I plotted this data I get the below map. enter image description here

However I wanted to convert these coordinates to lat/long and reclassify the Aragonite values into categories of 1 - 5. So I used the rasterToPoints(r, spatial = TRUE) function to create a SpatialPointsDataFrame that I reprojected to "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" using the spTransform function.

head(SPDF2)
  #Aragonite     long      lat
  #1 1.542071 20.89833 84.66917
  #2 1.538187 22.69496 84.66917
  #3 1.537830 24.49159 84.66917
  #4 1.534834 26.28822 84.66917
  #5 1.534595 28.08485 84.66917
  #6 1.532505 29.88148 84.66917

At this point I reclassified the data onto my 1-5 scale by writing a CSV file, reclassifying the Aragonite values and reading back into R.

head(r.csv)
      long      lat Aragonite
1 20.89833 84.66917        5
2 22.69496 84.66917        5
3 24.49159 84.66917        5
4 26.28822 84.66917        5
5 28.08485 84.66917        5
6 29.88148 84.66917        5

I have created a SpatialPointsDataFrame

coordinates(r.csv) <- ~ long + lat

# And assigned the correct CRS 
proj4string(r.csv) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#I then defined the vector I wanted to plot 
plotr <- r.csv$Aragonite

# Defined number of colours to be used in plot
nclr <- 5 #the number of categories for Aragonite
mypalette <- brewer.pal(nclr, "Reds")

# Plot the map with desired lat/long coordinates and data points with colour coding
map("world", fill=TRUE, col="white", bg="lightblue", ylim=c(-60, 90), mar=c(0,0,0,0))
points(r.csv$long, r.csv$lat, pch = 15, col= mypalette, cex = 0.5)

HOWEVER, the resulting map is a mess! The points seem to be in the right place but the colours don't seem to match the new Aragonite categories at all? enter image description here

I tried another tack and tried to rasterize my points but had no luck as they are irregularly spaced.

3

Why exactly do you call rasterToPoints in the first place? Raster-based reclassifications can easily be done using reclassify from the raster package. Here is a brief example using SST data from the remote package:

library(remote)
data("pacificSST")

## custom reclassification matrix
rcl <- matrix(c(0, 27.5, 1, 
                27.5, 30, 2, 
                30, 31, 3), ncol = 3, byrow = TRUE)

## reclassify the first SST layer
rcl <- reclassify(pacificSST[[1]], rcl = rcl, right = FALSE)

## visualize
library(RColorBrewer)
clr <- rev(brewer.pal(3, "RdYlBu"))

library(rworldmap)
library(latticeExtra)
spplot(rcl, col.regions = clr, at = seq(.5, 3.5, 1), 
       scales = list(draw = TRUE)) + 
  layer(sp.polygons(countriesCoarse))

spplot

Categorial rasters can be visualized nicely using levelplot from rasterVis, so you should probably have a look at it.

  • The reason I did raster to points was to then use spTransform to reproject the coordinates to latitude and longitude. I reprojected the spatial points data frame because I couldn't reproject the raster itself (see here: gis.stackexchange.com/questions/220589/…) So now I have my data in a csv. file as above that I cannot seem to plot – Phoebe Stewart-Sinclair Dec 8 '16 at 9:29
  • Yeah this is what I have been trying to do but have come up with multiple error messages that I can't seem to rectify. gis.stackexchange.com/questions/220589/… – Phoebe Stewart-Sinclair Dec 8 '16 at 9:35
  • Ah, sorry @PhoebeStewart-Sinclair. I didn't realize these two questions were connected. Anyway, now that the projection issue is solved, image classification should be a lot easier! – fdetsch Dec 9 '16 at 7:08
  • Will use your reclassification script now instead of doing manually in a csv file, thanks a lot! – Phoebe Stewart-Sinclair Dec 9 '16 at 12:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.