Skip to content

Commit

Permalink
old
Browse files Browse the repository at this point in the history
  • Loading branch information
@rhijmans
rhijmans committed Apr 23, 2024
1 parent c48b9e0 commit d141089
Show file tree
Hide file tree
Showing 8 changed files with 124 additions and 21 deletions.
14 changes: 8 additions & 6 deletions R/distance.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,14 +150,16 @@ test.for.lonlat <- function(xy) {
}

setMethod("distance", signature(x="matrix", y="matrix"),
function(x, y, lonlat, pairwise=FALSE) {
function(x, y, lonlat, pairwise=FALSE, unit="m") {
if (missing(lonlat)) {
lonlat <- test.for.lonlat(x) & test.for.lonlat(y)
warn("distance", paste0("lonlat not set. Assuming lonlat=", lonlat))
}
stopifnot(ncol(x) == 2)
stopifnot(ncol(y) == 2)
v <- vect()
stopifnot(unit %in% c("m", "km"))
m <- ifelse(unit == "m", 1, 0.001)
d <- v@ptr$point_distance(x[,1], x[,2], y[,1], y[,2], pairwise[1], 1, lonlat)
messages(v)
if (pairwise) {
Expand All @@ -169,14 +171,14 @@ setMethod("distance", signature(x="matrix", y="matrix"),
)

setMethod("distance", signature(x="data.frame", y="data.frame"),
function(x, y, lonlat, pairwise=FALSE) {
function(x, y, lonlat, pairwise=FALSE, unit="m") {
distance(as.matrix(x), as.matrix(y), lonlat, pairwise=pairwise)
}
)


setMethod("distance", signature(x="matrix", y="missing"),
function(x, y, lonlat=NULL, sequential=FALSE, pairs=FALSE, symmetrical=TRUE) {
function(x, y, lonlat=NULL, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m") {

if (missing(lonlat)) {
lonlat <- test.for.lonlat(x) & test.for.lonlat(y)
Expand All @@ -186,13 +188,13 @@ setMethod("distance", signature(x="matrix", y="missing"),
crs <- ifelse(isTRUE(lonlat), "+proj=longlat +datum=WGS84",
"+proj=utm +zone=1 +datum=WGS84")
x <- vect(x, crs=crs)
distance(x, sequential=sequential, pairs=pairs, symmetrical=symmetrical)
distance(x, sequential=sequential, pairs=pairs, symmetrical=symmetrical, unit=unit)
}
)

setMethod("distance", signature(x="data.frame", y="missing"),
function(x, y, lonlat=NULL, sequential=FALSE, pairs=FALSE, symmetrical=TRUE) {
distance(as.matrix(x), lonlat=lonlat, sequential=sequential, pairs=pairs, symmetrical=symmetrical)
function(x, y, lonlat=NULL, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m") {
distance(as.matrix(x), lonlat=lonlat, sequential=sequential, pairs=pairs, symmetrical=symmetrical, unit=unit)
}
)

Expand Down
6 changes: 3 additions & 3 deletions R/spatvec.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,7 +200,7 @@ setMethod("fillHoles", signature(x="SpatVector"),
} else {
x@ptr <- x@ptr$remove_holes()
}
messages(x)
messages(x, "fillHoles")
}
)

Expand All @@ -222,7 +222,7 @@ setMethod("centroids", signature(x="SpatVector"),
} else {
x@ptr <- x@ptr$centroid(TRUE)
}
messages(x)
messages(x, "centroids")
}
)

Expand All @@ -231,7 +231,7 @@ setMethod("centroids", signature(x="SpatVector"),
setMethod("densify", signature(x="SpatVector"),
function(x, interval, equalize=TRUE, flat=FALSE) {
x@ptr <- x@ptr$densify(interval, equalize, flat)
messages(x)
messages(x, "densify")
}
)

Expand Down
1 change: 1 addition & 0 deletions src/RcppModule.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -962,6 +962,7 @@ RCPP_MODULE(spat){
.method("sieve", &SpatRaster::sieveFilter)
.method("view", &SpatRaster::viewshed)
.method("proximity", &SpatRaster::proximity)
.method("fillNA", &SpatRaster::fillNA)

.method("rectify", &SpatRaster::rectify)
.method("stretch", &SpatRaster::stretch)
Expand Down
20 changes: 10 additions & 10 deletions src/distRaster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -847,7 +847,7 @@ std::vector<double> SpatVector::pointdistance(const std::vector<double>& px, con
if (szp == szs) {
if (lonlat) {
for (size_t i = 0; i < szs; i++) {
d.push_back( distance_lonlat(px[i], py[i], sx[i], sy[i]) );
d.push_back( distance_lonlat(px[i], py[i], sx[i], sy[i]) * m);
}
} else { // not reached
for (size_t i = 0; i < szs; i++) {
Expand All @@ -857,7 +857,7 @@ std::vector<double> SpatVector::pointdistance(const std::vector<double>& px, con
} else if (szp == 1) { // to avoid recycling.
if (lonlat) {
for (size_t i = 0; i < szs; i++) {
d.push_back( distance_lonlat(px[0], py[0], sx[i], sy[i]));
d.push_back( distance_lonlat(px[0], py[0], sx[i], sy[i]) * m);
}
} else { // not reached
for (size_t i = 0; i < szs; i++) {
Expand All @@ -867,19 +867,19 @@ std::vector<double> SpatVector::pointdistance(const std::vector<double>& px, con
} else { // if (szs == 1) {
if (lonlat) {
for (size_t i = 0; i < szp; i++) {
d.push_back( distance_lonlat(px[i], py[i], sx[0], sy[0]));
d.push_back( distance_lonlat(px[i], py[i], sx[0], sy[0]) * m);
}
} else { // not reached
for (size_t i = 0; i < szp; i++) {
d.push_back( distance_plane(px[i], py[i], sx[0], sy[0]) * m );
d.push_back( distance_plane(px[i], py[i], sx[0], sy[0]) * m);
}
}
}
} else {
if (lonlat) {
for (size_t i=0; i<szp; i++) {
for (size_t j=0; j<szs; j++) {
d.push_back(distance_lonlat(px[i], py[i], sx[j], sy[j]));
d.push_back(distance_lonlat(px[i], py[i], sx[j], sy[j]) * m);
}
}
} else { // not reached
Expand Down Expand Up @@ -974,7 +974,7 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::stri

std::string distfun="";
d = geos_distance(x, pairwise, distfun);
if ((!lonlat) && (m != 1)) {
if (m != 1) {
for (double &i : d) i *= m;
}
return d;
Expand All @@ -993,7 +993,7 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::stri
if (s == sx) {
if (lonlat) {
for (size_t i = 0; i < s; i++) {
d[i] = distance_lonlat(p[0][i], p[1][i], px[0][i], px[1][i]);
d[i] = distance_lonlat(p[0][i], p[1][i], px[0][i], px[1][i]) * m;
}
} else { // not reached
for (size_t i = 0; i < s; i++) {
Expand All @@ -1003,7 +1003,7 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::stri
} else if (s == 1) { // to avoid recycling.
if (lonlat) {
for (size_t i = 0; i < sx; i++) {
d[i] = distance_lonlat(p[0][0], p[1][0], px[0][i], px[1][i]);
d[i] = distance_lonlat(p[0][0], p[1][0], px[0][i], px[1][i]) * m;
}
} else { // not reached
for (size_t i = 0; i < sx; i++) {
Expand All @@ -1013,7 +1013,7 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::stri
} else { // if (sx == 1) {
if (lonlat) {
for (size_t i = 0; i < s; i++) {
d[i] = distance_lonlat(p[0][i], p[1][i], px[0][0], px[1][0]);
d[i] = distance_lonlat(p[0][i], p[1][i], px[0][0], px[1][0]) * m;
}
} else { // not reached
for (size_t i = 0; i < s; i++) {
Expand All @@ -1026,7 +1026,7 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::stri
for (size_t i=0; i<s; i++) {
size_t k = i * sx;
for (size_t j=0; j<sx; j++) {
d[k+j] = distance_lonlat(p[0][i], p[1][i], px[0][j], px[1][j]);
d[k+j] = distance_lonlat(p[0][i], p[1][i], px[0][j], px[1][j]) * m;
}
}
} else { // not reached
Expand Down
91 changes: 91 additions & 0 deletions src/gdal_algs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1973,6 +1973,97 @@ SpatRaster SpatRaster::fillna(int threshold, int connections, SpatOptions &opt)
CPLErr GDALFillNodata(GDALRasterBandH hTargetBand, GDALRasterBandH hMaskBand, doubledfMaxSearchDist, intbDeprecatedOption, intnSmoothingIterations, char**papszOptions, GDALProgressFuncpfnProgress, void*pProgressArg)
*/


SpatRaster SpatRaster::fillNA(double missing, double maxdist, int niter, SpatOptions &opt) {

SpatRaster out = geometry(1, true, true, true);

if (!hasValues()) {
out.setError("input raster has no values");
return out;
}
if (maxdist <= 0) {
out.setError("maxdist should be > 0");
return out;
}
if (niter < 0) {
out.setError("niter should be >= 0");
return out;
}

std::string filename = opt.get_filename();
std::string driver;
if (filename.empty()) {
if (canProcessInMemory(opt)) {
driver = "MEM";
} else {
filename = tempFile(opt.get_tempdir(), opt.tmpfile, ".tif");
opt.set_filenames({filename});
driver = "GTiff";
}
} else {
driver = opt.get_filetype();
getGDALdriver(filename, driver);
if (driver.empty()) {
setError("cannot guess file type from filename");
return out;
}
std::string errmsg;
if (!can_write({filename}, filenames(), opt.get_overwrite(), errmsg)) {
out.setError(errmsg);
return out;
}
}

SpatOptions ops(opt);
GDALDatasetH hSrcDS, hDstDS;
if (!open_gdal(hSrcDS, 0, false, ops)) {
out.setError("cannot open input dataset");
return out;
}

GDALDriverH hDriver = GDALGetDriverByName( driver.c_str() );
if ( hDriver == NULL ) {
out.setError("empty driver");
return out;
}

//opt.datatype = "INT4S";
if (!out.create_gdalDS(hDstDS, filename, driver, true, 0, source[0].has_scale_offset, source[0].scale, source[0].offset, opt)) {
out.setError("cannot create new dataset");
GDALClose(hSrcDS);
return out;
}

GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDS, 1);
GDALRasterBandH hTargetBand = GDALGetRasterBand(hDstDS, 1);

if (GDALFillNodata(hTargetBand, hSrcBand, maxdist, 0, niter, NULL, NULL, NULL) != CE_None) {
out.setError("fillNA failed");
GDALClose(hSrcDS);
GDALClose(hDstDS);
return out;
}

GDALClose(hSrcDS);
if (driver == "MEM") {
if (!out.from_gdalMEM(hDstDS, false, true)) {
out.setError("conversion failed (mem)");
}
GDALClose(hDstDS);
return out;
}

double adfMinMax[2];
GDALComputeRasterMinMax(hTargetBand, true, adfMinMax);
GDALSetRasterStatistics(hTargetBand, adfMinMax[0], adfMinMax[1], -9999, -9999);
GDALClose(hDstDS);
return SpatRaster(filename, {-1}, {""}, {}, {});
}




/*
#include <gdalpansharpen.h>
SpatRaster SpatRaster::panSharpen(SpatRaster pan, SpatOptions &opt) {
Expand Down
5 changes: 3 additions & 2 deletions src/geos_methods.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1154,7 +1154,8 @@ SpatVector SpatVector::buffer(std::vector<double> d, unsigned quadsegs, std::str
}
b[i] = geos_ptr(pt, hGEOSCtxt);
}
SpatVectorCollection coll = coll_from_geos(b, hGEOSCtxt);
const std::vector<long> ids = std::vector<long>();
SpatVectorCollection coll = coll_from_geos(b, hGEOSCtxt, ids, false);

GEOSBufferParams_destroy_r(hGEOSCtxt, bufparms);
geos_finish(hGEOSCtxt);
Expand All @@ -1165,6 +1166,7 @@ SpatVector SpatVector::buffer(std::vector<double> d, unsigned quadsegs, std::str
}



// basic version of buffer, for debugging
SpatVector SpatVector::buffer2(std::vector<double> d, unsigned quadsegs) {

Expand Down Expand Up @@ -2985,7 +2987,6 @@ SpatVector SpatVector::cross_dateline(bool &fixed) {




SpatVector SpatVector::centroid(bool check_lonlat) {

SpatVector out;
Expand Down
1 change: 1 addition & 0 deletions src/spatRaster.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -643,6 +643,7 @@ class SpatRaster {
SpatRaster disaggregate(std::vector<unsigned> fact, SpatOptions &opt);
SpatRaster distance(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, bool haversine, SpatOptions &opt);
SpatRaster proximity(double target, double exclude, bool keepNA, std::string unit, bool buffer, double maxdist, bool remove_zero, SpatOptions &opt);
SpatRaster fillNA(double missing, double maxdist, int niter, SpatOptions &opt);

SpatRaster distance_rasterize(SpatVector p, double target, double exclude, std::string unit, bool haversine, SpatOptions &opt);
SpatRaster direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, SpatOptions &opt);
Expand Down
7 changes: 7 additions & 0 deletions src/write_gdal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -570,6 +570,13 @@ bool SpatRaster::writeStartGDAL(SpatOptions &opt, const std::vector<std::string>
std::vector<std::string> nms = getNames();
double naflag=NAN;
bool hasNAflag = opt.has_NAflag(naflag);

// if (driver == "AAIGrid" && std::isnan(naflag)) {
// avoid nan as flag
// naflag = -3.40282347E+38;
// hasNAflag = true;
// set opt.NAflag?
// }

if (writeRGB) nms = {"red", "green", "blue"};

Expand Down

0 comments on commit d141089

Please sign in to comment.