user interface for contouring; r-spatial/stars#13

DESCRIPTION

@@ -47,7 +47,7 @@ Depends:
     methods,
     R (>= 3.3.0)
 Imports:
-    classInt,
+    classInt (>= 0.2-1),
     DBI (>= 0.8),
     graphics,
     grDevices,

NEWS.md

@@ -1,12 +1,14 @@
 # version 0.7-0
 
-* `read_sf` no longer first creates tibbles from `data.frame`s, but creates them directly; #853
+* `st_read` receives a `query` argument that can run queries against OGR datasets; #834, by Barry Rowlingson and Michael Sumner
+
+* `read_sf` no longer first creates tibbles from `data.frame`s, but creates them directly; #853, db propagation by Etienne Racine
 
 * check difference between compile-time and run-time GEOS versions; #844
 
-* all GEOS routines are now more robust against memory leaks, by using unique pointers; #822, #845, by Dan Baston
+* all GEOS routines are now (more) robust against memory leaks, by using unique pointers; #822, #845, by Dan Baston
 
-* `st_buffer` receives the buffer styles `endCapStyle`, `joinStyle` and `mitreLimit`; #833, #842 by Mike Sumner
+* `st_buffer` receives the buffer styles `endCapStyle`, `joinStyle` and `mitreLimit`; #833, #842 by Michael Sumner
 
 # version 0.6-4
 
@@ -447,7 +449,7 @@
 * rename `st_makegrid` to `st_make_grid`, and `st_linemerge` to `st_line_merge`
 * add NEWS.md file (#207)
 
-* faster conversion of `data.frame` into `POINT` `sf` object, using `st_as_sf` (Mike Sumner)
+* faster conversion of `data.frame` into `POINT` `sf` object, using `st_as_sf` (Michael Sumner)
 
 * `rbind` method for `sf` objects now keeps coordinate reference system
 

R/plot.R

@@ -132,12 +132,9 @@ plot.sf <- function(x, y, ..., main, pal = NULL, nbreaks = 10, breaks = "pretty"
 						values = as.numeric(values) # drop units, if any
 						if (is.character(breaks)) { # compute breaks from values:
 							n.unq = length(unique(na.omit(values)))
-							breaks = if (! all(is.na(values)) && n.unq > 1) {
-								if (utils::packageVersion("classInt") > "0.2-1")
-									classInt::classIntervals(na.omit(values), min(nbreaks, n.unq), breaks, warnSmallN = FALSE)$brks
-								else
-									classInt::classIntervals(na.omit(values), min(nbreaks, n.unq), breaks)$brks
-							} else
+							breaks = if (! all(is.na(values)) && n.unq > 1)
+								classInt::classIntervals(na.omit(values), min(nbreaks, n.unq), breaks, warnSmallN = FALSE)$brks
+							else
 								range(values, na.rm = TRUE) # lowest and highest!
 						}
 						# this is necessary if breaks were specified either as character or as numeric

R/stars.R

@@ -180,10 +180,13 @@ gdal_subdatasets = function(file, options = character(0), name = TRUE) {
 
 #' @param use_integer boolean; if \code{TRUE}, raster values are read as (and rounded to) unsigned 32-bit integers values; if \code{FALSE} they are read as 32-bit floating points numbers. The former is supposedly faster.
 #' @param mask stars object with NA mask (0 where NA), or NULL
+#' @param cuts cut values for contour polygons (or lines)
+#' @param use_contours logical;
+#' @param contour_lines logical;
 #' @name gdal
 #' @export
 gdal_polygonize = function(x, mask = NULL, file = tempfile(), driver = "GTiff", use_integer = TRUE,
-		geotransform, use_contours = FALSE, contour_options = character(0)) {
+		geotransform, cuts = classInt::classIntervals(x[[1]])$brks, use_contours = FALSE, contour_lines = FALSE) {
 	gdal_write(x, file = file, driver = driver, geotransform = geotransform)
 	on.exit(unlink(file))
 	mask_name = if (!is.null(mask)) {
@@ -193,9 +196,18 @@ gdal_polygonize = function(x, mask = NULL, file = tempfile(), driver = "GTiff",
 			mask_name
 		} else
 			character(0)
+	contour_options = if (use_contours) # construct contour_options:
+			c(paste0("FIXED_LEVELS=", paste0(cuts, collapse = ",")),
+			paste0("ELEV_FIELD=0"),
+			paste0("POLYGONIZE=", ifelse(contour_lines, "NO", "YES")))
+		else
+			character(0)
 	pol = CPL_polygonize(file, mask_name, "GTiff", "Memory", "foo", character(0), 0, contour_options, use_contours, use_integer)
 	out = process_cpl_read_ogr(pol, quiet = TRUE)
 	names(out)[1] = names(x)[1]
+	if (use_contours)
+		out[[1]] = structure(match(out[[1]], cuts) - 1, 
+			levels = levels(cut(cuts, cuts, include.lowest=TRUE)), class = "factor")
 	out
 }
 

src/polygonize.cpp

@@ -84,42 +84,40 @@ Rcpp::List CPL_polygonize(Rcpp::CharacterVector raster, Rcpp::CharacterVector ma
 	OGRLayer *poLayer = poDS->CreateLayer("raster", sr, wkbMultiPolygon, NULL);
 	delete sr;
 
-	// create field:
-	OGRFieldDefn oField("Name", OFTInteger);
-	if (poLayer->CreateField(&oField) != OGRERR_NONE)
-    	Rcpp::stop("Creating attribute field failed.\n");
-
 	CPLErr err;
-	if (use_contours) {
-#if ((GDAL_VERSION_MAJOR > 3) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 4)) 
-		OGRFieldDefn idField("ID", OFTInteger);
-		if (poLayer->CreateField(&idField) != OGRERR_NONE)
-    		Rcpp::stop("Creating attribute field failed.\n");
-		OGRFieldDefn elevField("ELEV", OFTReal);
-		if (poLayer->CreateField(&elevField) != OGRERR_NONE)
+	if (use_integer) {
+		// create field:
+		OGRFieldDefn oField("Value", OFTInteger);
+		if (poLayer->CreateField(&oField) != OGRERR_NONE)
     		Rcpp::stop("Creating attribute field failed.\n");
-		err = GDALContourGenerateEx((GDALRasterBandH) poBand, (void *) poLayer,
-                       create_options(contour_options).data(), NULL, NULL);
-		if (err != OGRERR_NONE)
-			Rcpp::Rcout << "GDALContourGenerateEx returned an error" << std::endl;
-#else
-		Rcpp::stop("contour only available in GDAL >= 2.4.0");
-#endif
-	} else if (use_integer)
-		err = GDALPolygonize((GDALRasterBandH) poBand, maskBand,
+		if (GDALPolygonize((GDALRasterBandH) poBand, maskBand,
 			(OGRLayerH) poLayer,
 			iPixValField[0],
 			NULL, // create_options(options, true),
-			NULL, NULL);
-	else
-		err = GDALFPolygonize((GDALRasterBandH) poBand, maskBand,
-			(OGRLayerH) poLayer,
-			iPixValField[0],
-			NULL, // create_options(options, true),
-			NULL, NULL);
+			NULL, NULL) != OGRERR_NONE)
+				Rcpp::Rcout << "GDALPolygonize returned an error" << std::endl;
+	} else {
+		OGRFieldDefn oField("Value", OFTReal);
+		if (poLayer->CreateField(&oField) != OGRERR_NONE)
+    		Rcpp::stop("Creating attribute field failed.\n");
 
-	if (err != OGRERR_NONE)
-		Rcpp::Rcout << "GDAL[F}Polygonize returned an error" << std::endl;
+		if (!use_contours) {
+			if (GDALFPolygonize((GDALRasterBandH) poBand, maskBand,
+				(OGRLayerH) poLayer,
+				iPixValField[0],
+				NULL, // create_options(options, true),
+				NULL, NULL) != OGRERR_NONE)
+					Rcpp::Rcout << "GDALFPolygonize returned an error" << std::endl;
+		} else {
+#if ((GDAL_VERSION_MAJOR > 3) || (GDAL_VERSION_MAJOR == 2 && GDAL_VERSION_MINOR >= 4)) 
+			if (GDALContourGenerateEx((GDALRasterBandH) poBand, (void *) poLayer,
+                       	create_options(contour_options).data(), NULL, NULL) != OGRERR_NONE)
+				Rcpp::stop("GDALContourGenerateEx returned an error");
+#else
+			Rcpp::stop("contour only available in GDAL >= 2.4.0");
+#endif
+		} 
+	}
 
 	Rcpp::NumericVector type(1);
 	type[0] = 0;