Implementation of Smart Sampling

The vagaries of weather, especially in the background of climate change, has forced Indian Government to take firm steps to relieve farmers from the related stress. The situation gets complicated when farmers are forced to commit suicide in the want of money due to failed crops. By keeping this in the view, the Indian Government brought a very ambitious insurance scheme namely: Pradhan Mantri Fasal Bima Yojana (PMFBY). This is not only an insurance scheme but it does push for adoption of new and cutting edge technologies such as Geographic Information System (GIS), Remote Sensing, Drones, Artificial Intelligence (AI) etc. for the better assessment of nation-wide agriculture conditions and productivity.

    Furthermore, PMFBY is a yield-based crop insurance scheme that uses Crop Cutting Experiments (CCEs) to determine the yield loss suffered by the farmers due to natural catastrophes and adverse weather conditions. As an improvement on the CCEs, various agencies are trying to perfect the method of yield calculation and reporting. One such method is to amalgamate CCE with technology to select the sample sites for CCEs. These samples must be representative of the crop conditions in a given area, normally a Gram Panchayat(GP). A method in this regarding was proposed by Murthy et al. (1997) whereby Normalized Differential Vegetation Index (NDVI) from satellite images is to be used for assessing the conditions (thus creating strata) of crop and then taking a random sample from each such strata. In a sense, smart sampling is nothing but stratified sampling method. Now, the question comes how to implement this. Since there are no off-the-shelf tools available for smart sampling, a tool need to be produced by combining other smaller tools. This can be done in ArcGIS, R, Python etc. platforms.

Here I am showing how smart sampling can be done in freely available statistical programming language R:

R code for Smart Sampling using NDVI

#First Libraries
library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.5-17, (SVN revision 1070)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/Gulshan/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/Gulshan/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

# read the NDVI data
ndvi <- raster("path_for_NDVI_image.tif")
gp <- readOGR("path_of_directory_whereyour_gp_locationis_stored", "shapefilename")

## OGR data source with driver: ESRI Shapefile 
## Source: "path shapefile will appear here", layer: "name of shapefile will appear here"
## with ## features
## It has ## fields
## Integer64 fields read as strings:  Tabulat_sh Tabulat__1

#If you NDVI image doesnot have the same projection as your shapefile, make them into one projection system.
# Reproject the NDVI image - one time process.

#newproj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#ndvi_repro <- projectRaster(ndvi, crs = newproj)

# Once the data is prepared loop the data through a for loop as below and it should work.



for(i in 1:dim(gp)[1]){
  a <- crop(ndvi_repro, extent(gp[i,]))
  b <- rasterize(gp[i,], a)
  a <- mask(x = a, mask = b)
  #here I am reclassifying the NDVI.
  class_diff <- (max(getValues(a), na.rm = TRUE) - min(getValues(a), na.rm = TRUE))/5
  z <- c(min(getValues(a), na.rm = TRUE), (min(getValues(a), na.rm = TRUE) + class_diff), 1,
      (min(getValues(a), na.rm = TRUE) + class_diff), (min(getValues(a), na.rm = TRUE) + 2 * class_diff), 2,
      (min(getValues(a), na.rm = TRUE) + 2 * class_diff), (min(getValues(a), na.rm = TRUE) + 3 * class_diff), 3,
      (min(getValues(a), na.rm = TRUE) + 3 * class_diff), (min(getValues(a), na.rm = TRUE) + 4 * class_diff), 4,
      (min(getValues(a), na.rm = TRUE) + 4 * class_diff), (max(getValues(a), na.rm = TRUE)), 5)
  reclass_z <- matrix(z, ncol = 3, byrow = TRUE)
  a <- reclassify(a, reclass_z, right = NA)
  x11() # this will pop up a new plotting window.
  plot(a)
  # if you wish to save the classified NDVI image; just remove # and give appropriate path.
  #writeRaster(a, filename = "givepathnameofclassifiedndvi", overwrite = TRUE)
  sampleSp <- sampleStratified(x = a, size = 1, xy = TRUE, sp = TRUE, na.rm=TRUE)
  sampleSp@data$area_per <- round((sampleSp@data$cell/sum(sampleSp@data$cell)) * 100, 1)
  plot(sampleSp, add = TRUE, pch = 18, col = 'red')
  title(i)
  # if you wish to save the ss points; just remove # and give appropriate path.
  #writeOGR(obj = sampleSp, dsn = "givedirectorypath_for_points_to_be_stored", layer = nameofthesamplepoints[i], driver="ESRI Shapefile", overwrite_layer = TRUE)
}

The points (selected smartly ;-)) will appear in red colour. Your output should look like the view in the video.

References:

Murthy, C. S.; Thiruvengadachari, S.; Jonna, S.; Raju, P. V. Design of Crop Cutting Experiments with Satellite Data for Crop Yield Estimation in Irrigated Command Areas. Geocarto International, 1997, 12 (2), 5–11.