diff --git a/R/hdist.R b/R/hdist.R
index 2d12e3dae85c3fc9210e1c554bb87b62cc1be3bf..a63a320d79e081b3753a3ffc2e69493d64cb5f34 100644
--- a/R/hdist.R
+++ b/R/hdist.R
@@ -104,12 +104,13 @@ hdist.sfc <- function(x, y, method="rghosh", gres=5, ditself=FALSE, maxsample=2.
     gdist <- loop_gdist(xdisc=xdisc, ydisc=ydisc, proj=proj, intersect=FALSE, parallel=parallel[2], cores=cores[2])
     if(ditself){
       if(verbose) cat("Computing Ghosh distance within catchments\n")
+      na_val <- set_units(NA, units(gdist), mode="standard")
       if(!identical(x,y)){
         gdist <- cbind(gdist,sapply(xdisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, intersect=FALSE, rescale=FALSE, diag=TRUE, parallel=parallel[2], cores=cores[2])}))
-        gdist <- rbind(gdist,c(sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, intersect=FALSE, rescale=FALSE, diag=TRUE, parallel=parallel[2], cores=cores[2])}),NA))
+        gdist <- rbind(gdist,c(sapply(ydisc,FUN=function(x){call_gdist(pts1=x, pts2=x, proj=proj, intersect=FALSE, rescale=FALSE, diag=TRUE, parallel=parallel[2], cores=cores[2])}),na_val))
       }else{
         gdist <- cbind(gdist,diag(gdist))
-        gdist <- rbind(gdist,c(diag(gdist),NA))
+        gdist <- rbind(gdist,c(diag(gdist),na_val))
       }
     units(gdist) <- units(st_distance(xdisc[[1]][1],ydisc[[1]][1])) #units are lost by some operations
     }