Image analysis in R: matching the TIFF image with the corresponding point pattern match
so I'm using spatstat package in R, and analyzing a TIFF segmentation mask image as well as the corresponding planar point pattern (using only the X,Y coordinates of the segmented object centroids).
For clarity, the planar point pattern object contains X,Y coordinates (centroids) of the objects from the TIFF segmentation image. I am trying to plot the X,Y coordinates in a layered manner which shows the mask and the centroids.
However when I load the TIFF image that contains a visual image of the boundary points the window. The dimensions of the TIFF image is 2800X2800
vim
integervalued pixel image
2800 x 2800 pixel array (ny, nx)
enclosing rectangle: [0, 2800] x [0, 2800] micrometer
The corresponding point patterns range from (X: 2X700) and Y: (2X700).
types$RB
Planar point pattern: 835 points
window: rectangle = [2.1364, 699.1579] x [1.5143, 699.3214] units (one unit = 2 micrometer)
How do I align up these two coordinate spaces in R using 'spatmat' or another image analysis package?
How do I align the coordinates of the centroids over the segmentation mask?
I could try to rescale the types$RB object so that it has the same dimensions as the TIFF image, however this does not really look like a match.
Update: For this approach, I think it is required to create a shape file from the TIFF file. and then begin overlaying the centroids (X,Y) coordinates over the shape file derived from a TIFF plot.
See also questions close to this topic

Calculare length of segments within a transect
I have a transect data with latitude, longitude and substrate types. Below I provide a script that creates a hypothetical data with 3 substrate types along a straight transect starting at longitude 24.5 and ending at 23.2. Within this transect there are 3 substrate types (a,b and c), but substrate type "a" occurs 4 times and substrate type "b" twice. I would like to calculate the total length (meters) of each "a","b" and "c" substrate type segments in the transect. As an example, the substrate segment "a" ends at the position of the first observation of "b" substrate type and the substrate segment c ends where the fourth "a" substrate type segment starts. I would like the length of. I have looked into the sp and Rdistance packages but I´m really stuck. With thanks in advance.
hypothetical example: each box denote each segment for which I would like to calculate the length of
Alon<c(23.20, 23.30,23.40,24.10,24.15, 23.95, 23.70, 23.60, 24.20, 24.25) Blon<c(23.80, 23.85, 24.00, 24.03, 24.06) Clon<c(23.47, 23.50,23.55) Alat<c(64,64,64,64,64, 64, 64, 64,64, 64) Blat<c(64,64, 64, 64,64) Clat<c(64,64, 64) A<as.data.frame(cbind(Alon, Alat)) B<as.data.frame(cbind(Blon, Blat)) C<as.data.frame(cbind(Clon, Clat)) plot(A$Alon, A$Alat, pch=97) points(B$Blon, B$Blat, col="red", pch=98) points(C$Clon, C$Clat, col="blue", pch=99) A$ID<seq.int(nrow(A)) A[,3]<"A" B$ID<seq.int(nrow(B)) B[,3]<"B" C$ID<seq.int(nrow(C)) C[,3]<"C" colnames(A)<c("lon", "lat", "ID") colnames(B)<c("lon", "lat", "ID") colnames(C)<c("lon", "lat", "ID") A<as.data.frame(A) B<as.data.frame(B) C<as.data.frame(C) pos< rbind(A,B,C) pos<pos[,c("ID","lon","lat")]

Preserve start points in UnionAggregate
DECLARE @Geom TABLE ( shape geometry, shapeType nvarchar(50) ); INSERT INTO @Geom(shape,shapeType) VALUES('LINESTRING(1 2, 3 4)', 'A'), ('LINESTRING(3.2 4, 7 8)', 'B'); SELECT * FROM @Geom SELECT geometry::UnionAggregate(shape).ToString(), geometry::UnionAggregate(shape) FROM @Geom;
The WKT for the output is
MULTILINESTRING ((7 8, 3.2 4), (3 4, 1 2))
when I would want
MULTILINESTRING ((1 2, 3 4), (3.2 4, 7 8))
Where the beginning of the "A" and "B" line should be
(1 2)
and(3.2 4)
respectfully.This behavior of
UnionAggregate
doesn't seem to care about "direction" of the geometry in order to maintain that A union B and B union A is the same result. However, I want to preserve start/endpoints as I am unioning street geometry and I want all the LINESTRINGs to go in their original direction.This problem is discussed here: https://social.msdn.microsoft.com/Forums/sqlserver/enUS/89e9536636494294a0bcf3921598157f/unionoflinestringsandreversingdirection?forum=sqlspatial
They seem to suggest at a possible solution about checking the end result, but it is not clear to me how to do that. It is hinted at in a linked thread that
The MultiLineString always represents the graph from the point which farthest from the origin point.
It is not clear to me exactly what this means, but I don't think I can just assume the result of a UnionAggregate is always the reverse of what I want
If it is hard to know directional intent then I can add M measures where the direction should follow increasing M values.
Assuming I have a method for reversing the points in line, how would I go about solving for this?
I found a function that mimics for
STUnion
for added support for Z and M measure: http://www.spatialdbadvisor.com/files/SQLServer.html#robo48 however it is noted that "their direction could change (eg Start/Start Point relationship).", which is what I want to avoid. 
Problem with correlog function in pgirness
I am trying to produce an autocorrelogram using output from the correlog function in pgirmess as is described here.
My UTM coordinate data are:
structure(list(V1 = c(698073L, 698095L, 698274L, 697806L, 698602L, 698632L, 697425L, 698272L, 698125L, 697935L, 698681L, 699287L, 698687L, 698042L, 698052L, 697477L, 698096L, 698782L, 698203L, 698113L, 698046L, 697923L, 699398L, 698143L, 698555L, 697973L, 698042L, 698080L, 697918L, 698253L, 698687L, 698719L, 698079L, 697982L, 698273L, 698995L, 698267L, 700678L, 698087L, 698887L, 698599L, 698883L, 697947L, 698159L, 697464L, 697897L, 698867L, 697775L, 698071L, 698043L, 698115L, 697979L, 697976L, 698064L, 698024L, 698048L, 698200L, 698092L, 698328L, 697998L, 697772L, 697957L, 698590L, 698789L, 698060L, 698071L, 699316L, 699094L, 698633L, 699289L, 698089L, 697819L, 697894L, 698000L, 700218L, 700004L, 699684L, 699743L, 726256L, 724864L, 725175L, 725438L, 725778L, 700386L, 700740L, 700701L, 724204L, 724049L, 723872L, 722994L, 698150L, 698197L, 698064L, 698285L, 698685L, 698716L, 698316L, 698170L, 698920L, 699335L, 698273L, 698100L, 698605L, 698089L, 697784L, 697693L, 697990L, 697809L, 697945L, 698108L, 697945L, 697736L, 698567L, 699035L, 697651L, 698062L, 699035L ), V2 = c(1853363L, 1853581L, 1853596L, 1854098L, 1853795L, 1853764L, 1853338L, 1853684L, 1853627L, 1853540L, 1853912L, 1853399L, 1853897L, 1853236L, 1853744L, 1853505L, 1853458L, 1853926L, 1853736L, 1853858L, 1853783L, 1853586L, 1853612L, 1853498L, 1853834L, 1853269L, 1853528L, 1853639L, 1853533L, 1853503L, 1853897L, 1853704L, 1853685L, 1853638L, 1853595L, 1853602L, 1853530L, 1852898L, 1853424L, 1853770L, 1853792L, 1853859L, 1853586L, 1853708L, 1853603L, 1852869L, 1853973L, 1853898L, 1853608L, 1853762L, 1853695L, 1853586L, 1853617L, 1853775L, 1853251L, 1853162L, 1853693L, 1853603L, 1853654L, 1853463L, 1853532L, 1853392L, 1853853L, 1853843L, 1853550L, 1853587L, 1853509L, 1853461L, 1853724L, 1853506L, 1853535L, 1853643L, 1853807L, 1853556L, 1852417L, 1852815L, 1853113L, 1852880L, 1785260L, 1786479L, 1786054L, 1785839L, 1785467L, 1852514L, 1852029L, 1852391L, 1787745L, 1787746L, 1787926L, 1788326L, 1853416L, 1853437L, 1853409L, 1853531L, 1853848L, 1853734L, 1853660L, 1853465L, 1854071L, 1853365L, 1853577L, 1853698L, 1853878L, 1853286L, 1853532L, 1853503L, 1853423L, 1854024L, 1854096L, 1853483L, 1853798L, 1853664L, 1853804L, 1852806L, 1853503L, 1853304L, 1852806L)), class = "data.frame", row.names = c(NA, 117L))
My response variable data are:
c(44.172, 43.975, 43.338, 43.854, 41.475, 42.414, 45.035, 44.025, 43.278, 44.116, 44.606, 44.242, 45.682, 44.495, 42.732, 44.02, 44.626, 43.535, 44.025, 43.409, 43.207, 43.46, 42.505, 45.348, 42.495, 44.505, 43.636, 42.879, 42.086, 43.364, 43.667, 43.005, 43.939, 44.813, 45.364, 42.475, 43.768, 42.909, 43.535, 44.949, 43.187, 42.5, 43.495, 43.318, 42.561, 43.747, 42.005, 44.293, 44.808, 43.702, 42.677, 44.566, 45.237, 43.859, 44.237, 43.909, 44.596, 44.667, 44.04, 44, 42.192, 44.929, 42.949, 43.944, 43.53, 43.298, 43.025, 43.424, 43.712, 42.51, 44.152, 44.01, 41.833, 43.505, 44.646, 44.566, 42.929, 43.636, 44.54, 43.232, 43.182, 42.217, 43.904, 42.864, 43.475, 42.323, 43.025, 43.429, 43.298, 44.914, 42.884, 43.601, 44.475, 42.505, 41.268, 44.227, 42.955, 42.833, 46.056, 44.293, 42.556, 43.717, 43.283, 43.086, 42.707, 43.859, 42.884, 44.682, 41.551, 42.803, 45.242, 42.929, 42.848, 44.288, 43.283, 43.288, 44.444)
I tried:
pgi.cor < correlog(coords=coords, z=sp.rich$mean], method="Moran",nbclass=11)
, but this appears to have cause problem calculating distances between the duplicate coordinates:Error in dnearneigh(coords[zero, ], breaks[i, 1], breaks[i, 2]) : nonpositive number of rows in x
Therefore, I tried:
pgi.cor < correlog(coords=unique(coords), z=sp.rich$mean[c(31,117)], method="Moran",nbclass=11)
(because rows 31 and 117 are duplicates of other rows). However this gives me the same error result. Any ideas? Thanks! 
Mannkendall trend analysis for Rainfall raster in R and its interpretation
I am doing a trend analysis for the raster. I need this plot.
The plot is MannKendalltrend analysis of annual precipitation... Colours show Sen's slope estimator units (mm y−1) for precipitation for each grid cell (1 km regular grid). Signiﬁcant trends (α= 0.1)are shown as hatched areas So far I have the following code.
##libraries library(raster) library(grDevices) library(spatialEco) ##download the data for tile with lat long datasets<getData('worldclim', var='prec', res=10) ##cropping to some random smaller extent e < as(extent(55,65,30,35), 'SpatialPolygons') crs(e) < "+proj=longlat +datum=WGS84 +no_defs" datasets < crop(datasets, e) names(datasets)<month.abb ##plot the data to check plot(datasets[[1]]) result<raster.kendall(datasets, tau = TRUE, intercept = TRUE, p.value = TRUE, z.value = TRUE, confidence = F) names(result) < c("slope","tau", "intercept", "p.value","z.value") ##make the significant p value shapefile p.value < result[[4]] p.value[p.value >= 0.05] < NA ##95%significance p.value<rasterToPolygons(p.value) p.value < aggregate(p.value) ##colour scale Uniques < cellStats(result[["slope"]],stat=unique) Uniques.max < max(Uniques) Uniques.min < min(Uniques) my.at < round(seq(ceiling(Uniques.max), floor(Uniques.min), length.out= 10),0) myColorkey < list(at=my.at, labels=list(at=my.at)) levelplot(result[["slope"]],margin=F,main="slope",par.settings=RdBuTheme(),colorkey=myColorkey,xlab=NULL,ylab=NULL,scales = list(draw=F,tck = 0))+ layer(sp.polygons(p.value,density=5,angle=45))
It gives this.
I need to verify the following.
Is this the correct way of computing the trend?
Can I interpret this slope as trend? Can I say the trend only for polygons at lower left corner are significant (<= 0.05). Even though the higher values of slope are at top right corner.
 How to fill the polygon as inclined hatches? I used density=5,angle=45 but it is not plotted.

How to migrate SP returning result sets from MySQL to Postgresql in simple way?
Currently I am having my business logic written in MySQL SP. Planning to migrate to Postgresql 11 (due to feature of multi threading in PG 11 unlike MySQL). There are nearly 300 SP consisting of SP with :
 SP Containing Transactions (simple update,insert and nothing to return back from SP).
 SP Returning multiple result sets (for instance 2 simple select queries). Whats the efficient way to migrate SP from MySQL to PG 11?
Solution: 1. We have tried PG 11 dynamic cursor in Functions to return multiple result sets. 2. Created simple SP (without returning result sets just update and insert) in PG till now.

Shiny and Leaflet integration is really slow  how can I speed it up?
Right now i'm almost certain that my current use of shiny and leaflet is suboptimal.
At a high level my current approach looks like this:
 Generate a leaflet.
 Create a reactive dataframe on user input.
 Create a reactive dataframe of lat lon coordinates on user selection of their area of interest.
 Merge a spatial dataframe (containing postcode polygon boundaries) with the reactive dataframe from step 2, then draw the map with the joined dataframe. This keeps all the data necessary for drawing polygons, adding colorBins and fillColor and labels inside the same final dataframe.
In more detail, the steps are executed as follows:
Generate a map like this:
output$leaflet_map < renderLeaflet({ leaflet() %>% addTiles() })
Produce a reactive dataframe of marketing data to be joined onto an
sf
spatial dataframe containing postcode polygons viasp::merge()
(the join happens a little later, i'll get to that):reactive_map_data1 < reactive({ df %>% filter(BrandAccount_Brand %in% input$selectBrandRecruitment1) %>% group_by(POA_CODE, ordertype) %>% summarise("Number of Orders type and postcode" = n(), "AOV" = round(mean(TotalDiscount), 2)) %>% left_join(seifa, by = "POA_CODE") %>% left_join(over25bypostcode, by = "POA_CODE") %>% mutate(`Proportion of Population Over 25` = round(n() / `25_and_over` * 100, 4)) })
Create a reactive dataframe containing the
lat
andlon
coordinates of the State selected by the user to be fed into the call to render the map:reactive_state_recruitment1 < reactive({ australian_states %>% filter(States == input$selectState_recruitment1) })
Render the final map 
profvis
determines that this is infact the slow part:observeEvent( input$gobutton_recruitment1, { ## First I load the spatial data with each call to render the ## map  this is almost certainly suboptimal however I can't ## think of another way to do this as each time the data are ## joined I have no other way of resetting the gdal.postcodes2 ## spatial dataframe to its original state which is why I reload ## it from .rds each time: gdal.postcodes_recruitment1 < readRDS("gdal.postcodes2.rds") ## I then merge the marketing `reactive_map_data1()` dataframe ## created in Step 2 with the freshly loaded `gdal.postcodes2` ## spatial dataframe  `profvis` says this is pretty slow but ## not as slow as the rendering of the map gdal.postcodes_recruitment1@data < sp::merge(gdal.postcodes_recruitment1@data, reactive_map_data1(), by.x = "POA_CODE", all.x = TRUE) ## Next I generate the domain of `colorBin` with the `Number of ## Orders type and postcode` variable that only exists after the ## merge and is subject to change from user input  it resides ## within the `reactive_map_data1()` dataframe that gets merged ## onto the `gdal.postcodes2()` spatial dataframe. pal < colorBin("YlOrRd", domain = gdal.postcodes_recruitment1$`Number of Orders type and postcode`, bins = bins_counts) ## Lastly I update the leaflet with `leafletProxy()` to draw the ## map with polygons and fill colour based on the ## `reactive_map_data1()` values leafletProxy("leaflet_map_recruitment1", data = gdal.postcodes_recruitment1) %>% addPolygons(data = gdal.postcodes_recruitment1, fillColor = ~pal(gdal.postcodes_recruitment1$`Number of Orders type and postcode`), weight = 1, opacity = 1, color = "white", dashArray = "2", fillOpacity = .32, highlight = highlightOptions( weight = 3.5, color = "white", dashArray = "4", fillOpacity = 0.35, bringToFront = TRUE), layerId = gdal.postcodes_recruitment1@data$POA_CODE, label = sprintf( "<strong>%s<br/>%s</strong><br/>%s<br/>%s<br/>%s<br/>%s", paste("Postcode: ", gdal.postcodes_recruitment1$POA_CODE, sep = ""), paste("% of Population Over 25: ", gdal.postcodes_recruitment1$`Proportion of Population Over 25`, "%"), paste("Number of Orders: ", gdal.postcodes_recruitment1$`Number of Orders type and postcode`, sep = ""), paste("Ave Order Value: $", gdal.postcodes_recruitment1$`AOV`, sep = ""), paste("Advantage & Disadvantage: ", gdal.postcodes_recruitment1$`Relative SocioEconomic Advantage and Disadvantage Decile`, sep = ""), paste("Education and Occupation: ", gdal.postcodes_recruitment1$`Education and Occupation Decile`, sep = "") ) %>% lapply(htmltools::HTML), labelOptions = labelOptions( style = list("fontweight" = "normal", padding = "3px 8px"), textsize = "15px", direction = "auto")) %>% addLegend("bottomright", pal = pal, values = ~bins_counts, title = "# of Recruits (All Time)", labFormat = labelFormat(suffix = ""), opacity = 1 ) %>% setView(lng = reactive_state_recruitment1()$Lon, lat = reactive_state_recruitment1()$Lat, zoom = reactive_state_recruitment1()$States_Zoom) })
All up the map takes between 7 and 20 seconds to render as the data are quite large.
Some things to note:
The polygons have already been simplified to death, they are currently only displaying at 10% of the detail that was originally provided to define postcode boundaries by the Australian Bureau of Statistics. Simplifying the polygons further is not an option.
sp::merge()
is not the fastest ofjoin
functions I have come across, but it is necessary in order to merge a spatial dataframe with a nonspatial dataframe (other joins such as those offered bydplyr
will not accomplish this task  a look at thesp::merge()
documentation reveals that this has something to do with S3 and S4 datatypes, in any case this part is not the slow part according toprofvis
).According to
profvis
the actual rendering of the map in step 4 (drawing polygons) is the slow part. Ideally a solution to speed this whole process up would involve drawing the polygons on the original leaflet, and only updating the fillColor and labels applied to each polygon upon input of the 'Go' actionButton. I have not figured out a way to do this.
Can anyone think of a way to restructure this whole procedure to optimise efficiency?
Any input is greatly appreciated.