Fix issue 11

DESCRIPTION

@@ -47,7 +47,6 @@ Suggests:
     pls,
     pscl,
     reshape2,
-    rgdal,
     rmarkdown,
     qicharts,
     qtl,

NEWS.md

@@ -4,7 +4,9 @@
 
 ## Other notes
 
-* Updated cochran.eelworms to fix a typo, and added more columns for grain yield, straw yield, weeds. Updates to documentation.
+* cochran.eelworms - Fix a typo reported by U.Genschel, and added more columns for grain yield, straw yield, weeds. Updates to documentation.
+
+* gartner.corn - Remove 'rgdal' package from example (Issue #11).
 
 # agridat 1.22 (2023-08-24)
 
@@ -20,7 +22,7 @@
 * hansen.multi.uniformity (multi-year uniformity trials)
 * heath.cabbage.uniformity
 * hutchinson.cotton.uniformity (Data from Rothamsted archive)
-* immer.sugarbeet.uniformity (add 1931 results from Rothamsted archive)
+* immer.sugarbeet.uniformity (add year 1931 results from Rothamsted archive)
 * khan.brassica.uniformity
 * kirk.potato (20 varieties with 15 reps)
 * mckinstry.cotton.uniformity (Data from Rothamsted archive)

man/gartner.corn.Rd

@@ -36,20 +36,21 @@
 
   Each harvested swath was 12 rows wide = 360 inches.
 
-  Time 0 is 5 Nov 2011, 12:38:03 Central Time.
-  Time 16359 = 4.54 hours later.
+  Timestamp 0 = 5 Nov 2011, 12:38:03 Central Time.
+  Timestamp 16359 = 4.54 hours later.
 
-  Yield is calculated as total dry weight (corrected to 15.5 percent moisture), divided by 56 pounds (to get bushels), divided by the harvested area.
+  Yield is calculated as total dry weight (corrected to 15.5 percent
+  moisture), divided by 56 pounds (to get bushels), divided by the
+  harvested area:
+
   drygrain = [massflow * seconds * (100-moisture) / (100-15.5)] / 56
   harvested area = (distance * swath width) / 6272640
   yield = drygrain / area
 
 }
 \source{
 
-  Originally from University of Minnesota Precision Agriculture Center.
-  https://www.soils.umn.edu/academics/classes/soil4111/hw/
-
+  University of Minnesota Precision Agriculture Center.
   Retrieved 27 Aug 2015 from https://web.archive.org/web/20100717003256/https://www.soils.umn.edu/academics/classes/soil4111/files/yield_a.xls
 
   Used via license: Creative Commons BY-SA 3.0.
@@ -61,7 +62,6 @@
   in on-farm experiments.
   Field Crops Research, 255, 15 September 2020, 107783.
   https://doi.org/10.1016/j.fcr.2020.107783
-
 }
 
 \examples{
@@ -71,21 +71,25 @@ library(agridat)
 data(gartner.corn)
 dat <- gartner.corn
 
-# Calculate yield
-dat <- transform(dat, yield=(mass*seconds*(100-moist)/(100-15.5)/56)/(dist*360/6272640))
-# Delete yield outliers
+# Calculate yield from mass & moisture
+dat <- transform(dat,
+  yield=(mass*seconds*(100-moist)/(100-15.5)/56)/(dist*360/6272640))
+
+# Delete low yield outliers
 dat <- subset(dat, yield >50)
 
-# Colors for yield
+# Group yield into 20 bins for red-gray-blue colors
 medy <- median(dat$yield)
 ncols <- 20
 wwidth <- 150
 brks <- seq(from = -wwidth/2, to=wwidth/2, length=ncols-1)
 brks <- c(-250, brks, 250) # 250 is safe..we cleaned data outside ?(50,450)?
 yldbrks <- brks + medy
 dat <- transform(dat, yldbin = as.numeric(cut(yield, breaks= yldbrks)))
+redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
+dat$yieldcolor = redblue(ncols)[dat$yldbin]
 
-# Add polygons for soil map units
+# Polygons for soil map units
 # Go to: https://websoilsurvey.nrcs.usda.gov/app/WebSoilSurvey.aspx
 # Click: Lat and Long. 43.924, -93.975
 # Click the little AOI rectangle icon.  Drag around the field
@@ -95,45 +99,44 @@ dat <- transform(dat, yldbin = as.numeric(cut(yield, breaks= yldbrks)))
 # Download the zip file and find the soilmu_a_aoi files.
 
 # Read shape files
-libs(rgdal)
-# shp <- readOGR("C:/x/rpack/agridat/inst/files","gartner.corn")
-shp <- readOGR(system.file(package="agridat", "files"), "gartner.corn")
-
-# Plot the shapefiles first to set up the coordinate system
-plot(shp, xlim=range(dat$long), ylim=range(dat$lat))
-box() # Add the yield points
-redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
-with(dat, points(long,lat, main="yield heat map",
-                 col=redblue(ncols)[yldbin], cex=.75, pch=16))
-plot(shp, add=TRUE, lwd=2) # Overlay soil polygons on top
-title("gartner.corn - yield heatmap with soil map unit symbols")
-# Manual annotation of soil map units
-text(x = c(-93.97641, -93.97787, -93.97550, -93.97693, -93.97654, -93.97480,
+libs(sf)
+fname <- system.file(package="agridat", "files", "gartner.corn.shp")
+shp <- sf::st_read( fname )
+
+# Annotate soil map units. Coordinates chosen by hand.
+mulabs = data.frame(
+ name=c("110","319","319","230","105C","110","211","110","211","230","105C"),
+ x = c(-93.97641, -93.97787, -93.97550, -93.97693, -93.97654, -93.97480,
        -93.97375, -93.978284, -93.977617, -93.976715, -93.975929),
-     y = c(43.92185, 43.92290, 43.92358, 43.92445, 43.92532, 43.92553,
-       43.92568, 43.922163, 43.926427, 43.926993, 43.926631),
-     lab=c("110","319","319","230","105C","110","211","110","211","230","105C"))
+ y = c(43.92185, 43.92290, 43.92358, 43.92445, 43.92532, 43.92553,
+       43.92568, 43.922163, 43.926427, 43.926993, 43.926631) )
+mulabs = st_as_sf( mulabs, coords=c("x","y"), crs=4326)
+mulabs = st_transform(mulabs, 2264)
 
-# Trim off the ends of the field & re-do image above
-dat <- subset(dat, lat < 43.925850 & lat > 43.921178)
+# Colored points for yield
+dat <- st_as_sf(dat, coords=c("long","lat"), crs=4326)
 
-# Identify the soil type for each yield point
-dat$ix <- over(SpatialPoints(dat[ , c('long','lat')]),
-               SpatialPolygons(shp@polygons))
-dat$mu <- shp@data[, "MUSYM"][dat$ix]
-# Check the points are properly identified
-# with(dat, points(long,lat, col=redblue(ncols)[ix], cex=.75, pch=16))
+libs(ggplot2)
+
+ggplot() +
+  geom_sf(data=dat, aes(col=yieldcolor) ) +
+  scale_color_identity() +
+  geom_sf_label(data=mulabs, aes(label=name), cex=2) +
+  geom_sf(data=shp["MUSYM"], fill="transparent") +
+  ggtitle("gartner.corn") +
+  theme_classic()
 
-# Aggregate points by soil type and analyze
-tapply(dat$yield, dat$mu, mean)
-tapply(dat$yield, dat$mu, sd)
-libs(lattice)
-densityplot(~yield|mu, dat, layout=c(1,5),
-            main="gartner.corn - yield density plot")
+
+# Further analyses...
+# Trim top and bottom ends of the field
+dat <- subset(dat, lat < 43.925850 & lat > 43.921178)
+# Identify the soil type for each yield point
+# Aggregate points by soil type
 
 
 if(0){
   # Draw a 3D surface.  Clearly shows the low drainage area
+  # Re-run the steps above up, stop before the "Colored points" line.
   libs(rgl)
   dat <- transform(dat, x=long-min(long), y=lat-min(lat), z=elev-min(elev))
   clear3d()