20 Nov

Hey,

Today a bit of trSIGonometry with this function to compute the angle/orientation/direction between 2 points of a GPX file.

```############################ #Function ############################   angle <- function(pt1, pt2) { conv = 180 / pi # conversion radians / degrees diff <- pt2 - pt1 diff.abs <- abs(diff) if (diff[1]>=0 & diff[2]>=0){angle <- pi/2 - atan(diff.abs[2] / diff.abs[1])} if (diff[1]>=0 & diff[2]<=0){angle <- pi/2 + atan(diff.abs[2] / diff.abs[1])} if (diff[1]<=0 & diff[2]>=0){angle <- 3*pi/2 + atan(diff.abs[2] / diff.abs[1])} if (diff[1]<=0 & diff[2]<=0){angle <- pi + (pi/2- atan(diff.abs[2] / diff.abs[1]))} return(angle * conv) }   ############################ #Script ############################   library(rgdal) library(raster) library(data.table)   #Set folder setwd("/home/GuruBlard/gpx/")   # Load file names filenames <- list.files(pattern=".gpx\$")   # Iteration through all gpx files for (filename in filenames){ # Load the GPX   #ogrListLayers("export.gpx") gpx.trackpoints <- readOGR(dsn = filename, layer= "track_points", stringsAsFactors = F)   # Extract the coordinates as a matrix: gpx.tp <- coordinates(gpx.trackpoints)   # #Loop through and add angle - the difference between the of the current point and the following one: for (tp in 1:(nrow(gpx.tp)-1)) { gpx.trackpoints\$angle[tp] <- angle(gpx.tp[tp,], gpx.tp[tp+1,]) }   gpx.trackpoints@data\$Date <- substr(as.character(gpx.trackpoints@data\$time), 1,10) gpx.trackpoints@data\$Hour <- substr(as.character(gpx.trackpoints@data\$time),12,19) gpx.trackpoints@data\$filename <- filename writeOGR(gpx.trackpoints, paste0(tools::file_path_sans_ext(filename), ".shp"), layer=tools::file_path_sans_ext(filename), driver="ESRI Shapefile") }```

17 Jul

Holà GISpo,

Today, a simple function to (highly) speed up the standard Euclidean distance function in R when you use large matrix (n> 100.000.000).

```euc_dist <- function(m) {mtm <- Matrix::tcrossprod(m); sq <- rowSums(m*m); sqrt(outer(sq,sq,"+") - 2*mtm)}   # Example M <- matrix(rexp(1000000, rate=.1), ncol=1000)   ptm <- proc.time() d <- dist(M, method = "euclidean") proc.time() - ptm   # user system elapsed # 1.08 0.00 1.08   ptm <- proc.time() d2 <- euc_dist(M) d2 <- as.dist(d2) proc.time() - ptm   # user system elapsed # 0.424 0.008 0.429   isTRUE(all.equal(as.matrix(d), as.matrix(d2))) #TRUE```
18 Apr

Hi PubliGISer,

You want to publish the revision of your awesome paper but you need a track change to meet reviewers’ requirements and you do not know how to get it? You are desperate and think that you will not be able to realize your dream of a successful research career?

Online LaTeX diff tool is the TOOL you are looking for.

Just copy paste the old and the new latex files and this tool will automatically compute the difference between the 2 documents.

Problem solved.

12 May

Hey folks ! It’s been a while !

Interested in meteorological data over US ? Here is a way to easily download Daymet data inspired from DaymetR package.

Daymet is a collection of algorithms and computer software designed to interpolate and extrapolate from daily meteorological observations to produce gridded estimates of daily weather parameters. Weather parameters generated include daily surfaces of minimum and maximum temperature, precipitation, humidity, and radiation produced on a 1 km x 1 km gridded surface over the conterminous United States, Mexico, and Southern Canada.

```setwd("/guru-blard/DAYMET")   param = c("vp", "tmin", "tmax", "swe", "srad", "prcp", "dayl") #see here https://daymet.ornl.gov/overview.html for the available variable tiles=c(11922:11925, 11742:11745) #id of the tiles of interest (cfr shapefile) year_range=1985:2015   for (i in year_range) { for (j in tiles) { for (k in param) { download_string = sprintf("https://daymet.ornl.gov/thredds/fileServer/ornldaac/1219/tiles/%s/%s_%s/%s.nc", i, j, i, k) daymet_file = paste(k, "_", i, "_", j, ".nc", sep = "") cat(paste("Downloading DAYMET data for tile: ", j, "; year: ", i, "; product: ", k, "\n", sep = "")) try(download(download_string, daymet_file, quiet = TRUE, mode = "wb"), silent = FALSE) } } }```
21 Jan

Happy New GIS,

Be sure to get the last version of GDAL (Ubuntu installation)

Go to Geofabrik.

Using ogr to convert the pbf file into a SQlite database which can be used directly in QGIS with Add SpatiaLite Layer.

`ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat xmin ymin xmax ymax GuruBlard_Region.db OSM_continent.pbf`

where `xmin ymin xmax ymax`:
spatial query extents, in the SRS of the source layer(s) (or the one specified with -spat_srs). Only features whose geometry intersects the extents will be selected. The geometries will not be clipped unless -clipsrc is specified

For further information, visit this post of Anita Graser.

31 Aug

Hello GIStouille,

Here is a nice function to merge two dataframe while keeping the original rows’ order (source).
keep_order=1 keeps order of the first dataframe (x).
keep_order=2 keeps order of the second dataframe (y).

```merge.with.order <- function(x,y, ..., sort = T, keep_order) { # this function works just like merge, only that it adds the option to return the merged data.frame ordered by x (1) or by y (2) add.id.column.to.data <- function(DATA) { data.frame(DATA, id... = seq_len(nrow(DATA))) } # add.id.column.to.data(data.frame(x = rnorm(5), x2 = rnorm(5))) order.by.id...and.remove.it <- function(DATA) { # gets in a data.frame with the "id..." column. Orders by it and returns it if(!any(colnames(DATA)=="id...")) stop("The function order.by.id...and.remove.it only works with data.frame objects which includes the 'id...' order column")   ss_r <- order(DATA\$id...) ss_c <- colnames(DATA) != "id..." DATA[ss_r, ss_c] }   # tmp <- function(x) x==1; 1 # why we must check what to do if it is missing or not... # tmp()   if(!missing(keep_order)) { if(keep_order == 1) return(order.by.id...and.remove.it(merge(x=add.id.column.to.data(x),y=y,..., sort = FALSE))) if(keep_order == 2) return(order.by.id...and.remove.it(merge(x=x,y=add.id.column.to.data(y),..., sort = FALSE))) # if you didn't get "return" by now - issue a warning. warning("The function merge.with.order only accepts NULL/1/2 values for the keep_order variable") } else {return(merge(x=x,y=y,..., sort = sort))} }```
25 Aug

Hi crazyGISfans,

Hdf format may be sometimes difficult to manage. We prefer working with geotif.
How do you do convert hdf in tif in R with gdal_translate?

```library(gdalUtils) library(tools)   setwd("/guru-blard/")   files<-list.files(, pattern=".hdf\$") band<-1   for (file in files){ gdal_translate(src_dataset=get_subdatasets(file)[band], dst_dataset=paste0(file_path_sans_ext(file),"_B", band,'.tif')) }```

That’s it.

08 Jun

Holà GISpa,

Here is a R function that uses the montage function of ImageBrick to assemble png in one grid. All the parameters are explained here.

```pattern='CoolestPngEver_' In<-list.files("/export/homes/gurublard/", full.names=T, pattern=pattern) Out<-paste0(/export/homes/gurublard/", pattern,"all.png") command<-paste0("montage -label '%f' -density 300 -tile 4x0 -geometry +5+50 -border 10 ", paste(In, collapse=" "), " ", Out) system(command)```

04 Mar

Dear followers,

Have you already need to list filenames in a directory that are different from a specific pattern ?

Just let the MaGIS happen with the list.diff function !

```list.diff<-function(dir, pattern){ filename.all<-list.files(dir) filename.diff<-list.files(dir, pattern) filename<-setdiff(filename.all,filename.diff) return(filename) }```