Land_cover_type_plot

Plot Land Cover Type

This document summaries how to plot a land cover map from a NetCDF file

# import NetCDF file as raster format
inputfile <- c("I:/Rscripts/migration_map/1994.nc")
# load library 
library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

library(raster)
library(rgdal)

## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/ychen/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/ychen/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-3

# Grab the lat and lon from the data
lat <- raster(inputfile, varname="nav_lat")

## Loading required namespace: ncdf4

lon <- raster(inputfile, varname="nav_lon")

# Convert to points and match the lat and lons
offset <- 0.008
plat <- rasterToPoints(lat)+ offset
plon <- rasterToPoints(lon)- offset
lonlat <- cbind(plon[,3], plat[,3])
# Specify the lonlat as spatial points with projection as long/lat
plonlat <- SpatialPoints(lonlat, proj4string = CRS("+proj=longlat +datum=WGS84"))

# show the plot 
ld_go <- TRUE
if (ld_go) {
  lu_plot <- as.matrix(raster(inputfile, varname="LU_TYPE"))
  LU_P <- lu_plot
  LU_P[LU_P==0 | LU_P >16  ] <- 0
  LU_P[LU_P==1]  <- 1.2
  LU_P[LU_P==4]  <- 2.2 
  LU_P[LU_P==2]  <- 3.2
  LU_P[LU_P==13] <- 4.2  
  LU_P[LU_P==8]  <- 5.5
  LU_P[LU_P>6 | LU_P<0] <-0
  
  breakpoints <- c(0,1,2,3,4,5,6)
  colors <- c("lightgray","#006400","#FFB90F","#A2CD5A","#0000FF","#FF1493")

  layout(matrix(data=c(1),nrow=1, ncol=1,  byrow=TRUE),
         widths=c(1,1,1), heights=c(1,1,1))
  
  par( oma=c(2,2,0,2), mgp=c(2,1,0), mar=c(3.0,3.0,0.0,2.0),plt = c(0.1,0.9,0.1,0.90))
  
  # assign the porjection 
  r_p <-raster(LU_P,
               xmn=extent(plonlat)[1], xmx=extent(plonlat)[2],
               ymn=extent(plonlat)[3],  ymx=extent(plonlat)[4], 
               crs=CRS("+proj=longlat +datum=WGS84") )
  plot(r_p, breaks=breakpoints,col=colors,
       lab.breaks=c("N/A","Forest","Agri.",
                    "Grass","to Water","B/Soil",""),
       xlab="",  ylab="", 
       xlim=c(119.8,122),ylim=c(21.8,25.3),cex=0.8)
  # add magrin text
  
  mtext("Longitude (deg)",side=1,line=2.5)
  mtext("Latitude (deg)",side=2,line=2.5)
  p_text <- paste("year 2015" ) 
  mtext(p_text,side=3,line=1.0,cex=2.0)
}