Tag Info

New answers tagged

1

I believe that you need to set the argument "longlat=TRUE". Since your data is in geographic coordinates it is likely that the kernel is being incorrectly defined. You also may want to explicitly specify the data slot "data = spdf@data". Please use caution with specification of the GWR method in anything other than exploratory analysis of nonstationarity. ...


2

You said: I am trying to extract values of raster cells on points using extract from raster package in R (similar to 'Extract Values to Points' in ArcGIS-10.2). Lets set up a test raster: require(raster) r = raster(ncol=100,nrow=100,xmn=0,xmx=1,ymn=0,ymx=1) r[]=runif(100*100) You said: After doing so in order to check integrity, I computed ...


-1

You might want to give a try to the GMT tool grdtrack. You can start by reconverting to your raster to a NetCDF file using gdal_translate: gdal_translate -of NetCDF myraster.tif myraster.grd and then put your x,y locations in a text file: x1, y1 x2, y2 .., .. then call grdtrack on that NetCDF grid based on your table of locations: grdtrack ...


0

I found an answer. I basically take everything I need into the colorkey function as follows. plot3 <- spplot(mymap, "percentage", col.regions = plotclr, at = round(class$brks, digits=1), colorkey = list(labels = list( labels = c("0%", "1%","2%","3%","4%","5%","6%","7%"), width = 2, cex = 2))) Width changes ...


0

You can define the labelling style using the label.style argument of lattice::contourplot. In my opinion, you should choose align. It is not a complete solution because it does not break the contour lines, but it is better than the default method. On the other hand you can overlay two different contour plots with different cuts, line widths, and labels ...


1

It's because NetCDF does not have names for each slice in the 3rd (and higher dims), but raster does. NetCDF has a name for a "variable" (which is the array), but raster has a name for every slice in the variable. (This is the standard mess where we conflate data fields/attributes with dimensions). There's no straightforward way to store these names in ...


3

This sorts itself out with the base function mean. mean(c(20,10),na.rm=TRUE) # where both values occur mean(c(20,NA),na.rm=TRUE) # where the first value occurs mean(c(NA,10),na.rm=TRUE) # where the second value occurs mean(c(NA,NA),na.rm=TRUE) # where both values are nodata If you think of raster functions in terms of vectorization then things become ...


3

As many on this forum know, I am often for an R solution. However, in this case it is reinventing the wheel, and in a much less robust way. There is a great piece of free software, Map Comparison Kit (MCK), that implements many published and novel validation statistics for rasters. Of particular interest in this case are the Kappa, fuzzy Kappa and weighted ...


0

If you are happy with the resolution provided by the GADM database, then you can easily download that data using the convenience function raster::getData(). (Do note: In addition to loading US county-level data as a SpatialPolygonsDataFrame to the current R session, the code below will also save the data for future use to a file named USA_adm2.Rdata and ...


0

This worked for me - the solution was to use SpatialPixelsDataFrame with the suggested tolerance argument (0.916421 in your case): points <- SpatialPoints(s100_ras[,c('x','y')], s100_ras[,c('z')] pixels <- SpatialPixelsDataFrame(points, tolerance = 0.916421, points@data) raster <- raster(pixels[,'z']) though, due to the high tolerance value, the ...


3

Why not just use county shapefiles from the Census Bureau? The Census recognizes those independent Virginia cities as both cities and as independent counties for the purposes of geographic distinction. Here is a list of the areas designated at the county level by the Census Bureau for the state of Virginia. You can find shapefiles for all counties here.


0

You can use spDistsN1 and which.min to assign the coordinates of the nearest points in another spatial points object. Things will remain ordered so, you could just assign the data from the original data back to the adjusted points. Add package and create offset points require(sp) data(meuse) pts <- meuse[1:10,] pts2 <- data.frame( ...


1

Here's a reproducible answer and a function that I think solves the problem. It all relies on nncross from the spatstat package. Step 1: Load the packages we'll be using library(sp) library(spatstat) library(maptools) # to convert to ppp Step 2: Create two small sets of points, give one attribute data: set.seed(2014) # ensure reproducibility x <- ...


4

Similar to @AndreJ, but use a dynamic gnomic projection, I mean a dynamic azimuthal equidistant projection for even more accuracy. An AEQ projection centred on each point will project equal distances in all directions, such as a buffered circle. (A Mercator projection will have some distortions in north and eastern directions, since it wraps around the side ...


4

You shall have a look at http://colorbrewer2.org/, and diverging color schemes therein. It can be imported into R using brewer.pal from the RColorBrewer package. library(RColorBrewer) cols <- brewer.pal(11, "RdBu") Or, for more than the default maximum of 11 colors, use colorRampPalette. cols <- colorRampPalette(brewer.pal(11, "RdBu"))


0

You can use the maptools::elide function to translate coordinates of a Spatial*object. If you combine the result with the sp.layout argument of spplot or with the latticeExtra::layer function, you will get what you need: library(sp) library(maptools) library(raster) library(lattice) library(latticeExtra) data(meuse) coordinates(meuse)=~x+y ZOOM <- ...


4

Instead of searching for the right UTM zone, you could create a custom transverse mercator projection for every point with +proj=tmerc +lat_0=.... +lon_0=... +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Draw the circle in that projection. The projected circle vertex coordinates will always be the same, so you have to create ...


0

What if you take the approach of creating a 1000 meter in EPSG:4326 around each of your points. Then convert the EPSG:4326 to your other coordinate system? The advantage of projecting the point, is that you don't have to worry about curvature of the earth with EPSG:4326.


3

I would write an R script that worked as a client, but will run on the database server. This will save the complication of trying to hook into PostGIS's backend and using PL/R (as I said in comments). The script will look something like this (which is practically pseudo-code here): > con = dbConnect(PG,"localhost","weather") # connect to local DB ...


1

the str( ) will display the internal structure of your object str(zambiap) If you just want a data.frame to export: dt<-data.frame(1:length(zambiap)) zambiapd <- SpatialPolygonsDataFrame(zambiap, data=dt) If you want a meaningful data.frame, you need to get the data from the zambiap object. dt.f<-NULL for (i in 1:length(zambiap)){ ...


1

For illustration sake, let's take a look at what rowSums, a special instance of apply, is actually doing. Understanding apply functions in R opens a door that allows you to run more complex analysis, optimize code and speed up processing. The reason that I expand on @cengal answer is because if one wants to write more complex functions that operate on rows ...


0

A spatial data frame works like a data frame. So if I understand your question correctly, you can just use rowSums. set.seed(107) mydf <- data.frame(var1 = c(sample(c(0:9), 10, replace = T)), var2 = c(sample(c(0:9), 10, replace = T)), var3 = c(sample(c(0:9), 10, replace = T)) mydf$NonZeroFields <- rowSums(mydf > 0) mydf var1 var2 var3 ...


1

Try: polygone1 <- gBuffer(polygone1, byid=TRUE, width=0) polygone2 <- gBuffer(polygone2, byid=TRUE, width=0) clip2 <- gIntersection(Polygone1, Polygone2, byid=TRUE) It is ugly, but it usually solves this kind of problem. HTH


2

Here is a suggestion using ggplot. I use ggplotGrob to combine the full and zoomed map and grid.arrange from the gridExtra add-on to combine the maps for different variables. There are many adjustments that can be made, of course. library(sp) library(ggplot2) library(grid) # for unit library(gridExtra) # for grid.arrange # zoom bounding box xlim <- ...


4

I'd say this isn't really a GIS question, and next time you might be better off posting it on Stack Overflow. Anyway, here is a quick example of how you might parse this using Python. import os, sys #directory holding the input files inDir = "C:\Temp" #Find all text files in the specified directory for inFile in os.listdir(inDir): if ...


6

Assuming all your .img files are in one folder and the ENVI file format suffices as the binary output this R code works. library(raster) directory <- "/path/to/IMG/files" setwd(directory) ## create vector containing all image filenames Images <- dir(directory, pattern="\\.img$") #edit the pattern (case sensitive!) if you want to exclude/include ...


2

You are actually asking two questions. One is how to merge census data to your shape file. This is not very difficult, you may want to check out merge or spCbind, or this answer. It hard to say without knowing the structure of your data. The second question is about how to do the weighted sampling. I have a suggestion below that uses spsample from the sp ...


1

You do not need a for loop. Just intersect everything at once and then add line lengths to the new line segments using the "SpatialLinesLengths" function in sp. Then, using the raster package rasterize function with the fun=sum argument you can create a raster with the sum of the line length(s) intersecting each cell. Using the above answer and associated ...


0

I guess that my answer from your other StackOverflow question did not lead you in the right direction? Here is a more detailed answer that may nudge you in the right direction. First, we need to know the projection of your data and the extent to be able to project to the Long/Lat grid correctly. Unfortunately, we do not have the PROJ4 CRS or the extent, so ...


1

That should be easy to do in QGIS. Load the postal Shapefile and join the sales rep data, then style based on sales rep. You can use the OpenLayers plugin to load background tiles (note though that these won't print).


1

Here is an approach in R using paste to concatenate fields and %in% to create an index of non-matching values that can be used to subset the data. You have some missing values in the "Feature_dataset_name" field so, I replaced them with "NA" values and this is propagated in the results (missing). setwd("D:/TEST") # Read data and concatenate fields f1 <- ...


0

If you need high-precision and robust geodesic measurements, use GeographicLib, which is natively written in several programming languages, including C++, Java, MATLAB, Python, etc. See C. F. F. Karney (2013) "Algorithms for geodesics" for a literary reference. Note that these algorithms are more robust and accurate than Vincenty's algorithm, for instance ...


2

This appears to have been a bug with the r 'raster' package version 2.2-31 released a bit before this post. I was able to work around the issue by upgrading the package to 2.3-12 which ran the exact same code without the error.


1

The best resource is the ggmap and ggplot2 packages in R. Here is a short paper describing details. http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf There are also materials from a workshop on visualizing climate data at this site: http://streaming.stat.iastate.edu/~dicook/NCAR/


0

Instead of points use SpatialPoints. Then you can reproject layers as you wish using spTransform from the rgdal package. library(maps) library(maptools) library(sp) library(rgdal) US <- map("state",fill=TRUE, plot=FALSE) US.names <- US$names US.IDs <- sapply(strsplit(US.names,":"),function(x) x[1]) US_poly_sp <- ...


0

Here is a suggestion using R. I have used the over function from the sp package to extract the points that fall within the polygons of the bounding boxes and assign the respective ID of the bounding box. library(sp) #read in points from csv input <- read.csv("input.csv", header=F) names(input) <- c("x","y","z") # create spatial points object pts ...


1

One possible way I can think of to accomplish this, using ArcGIS, is as follows: Add a field to each file, and calc its value to be a cancatenation of the other fields. Join file B to file A based on those cancatenated fields. Records from A which didn't have a record joined to it are unique to A (sort and select manually to get these, or do select by ...


0

With the updated question, the working directory has the name "LAS", which tells me this is LiDAR data. If you have LAS files, use these. Otherwise, convert your files to this much more efficient format using LAStools or libLAS (via OSGeo4W on Windows). First use txt2las on a headerless ASCII file, parsing X, Y, and Z fields: txt2las --verbose -parse xyz ...


0

You may want to explore the "rts" package. It specifically coerces raster stacks into time series objects. Although, I do not know if they are compatible with bfast. Your best bet would be to vectorize the problem in a function, passed to overlay, and use @Aaron advice by then coercing into a ts object. This is how I approach time series analysis of ...


5

Following a recent question, you may want to make use of the functionalities offered by the rgeos package to solve your problem. For reasons of reproducibility, I downloaded a shapefile of Tanzanian roads from DIVA-GIS and put it in my current working directory. For the upcoming tasks, you will need three packages: rgdal for general spatial data handling ...


6

If you are open to using alternative software to solve your problem, then I can suggest the Remove Off-Terrain Objects tool of the cross-platform open-source GIS Whitebox Geospatial Analysis Tools (of which I am lead developer). I realize that you said in your question that you could not convert your data to LAS format, but the tool takes a raster, not LAS ...


3

I would have to say that without the original LAS point cloud, you will only be incorporating more inaccuracies into the data through raster manipulation. The DEM provided looks relatively clean for a heavily urbanized 1m resolution DEM. The "uplifted squares" are a result of triangulations across data voids, where the buildings are not included in the ...


2

Not sure about the specific boundary that you are after but I believe that you are referring to the "maps" package. If there is not adequate detail in the maps datasets, there is a supplementary package "mapdata" with a worldHires dataset. From the package description: "This world database comes from a cleaned-up version of the CIA World Data Bank II data ...


5

I found a way using purely Python to get the coordinates for tweets using a word filter. It doesn't seem like many people include location with their tweets. This might not be what you're after either because this is live streaming data. You can test it by putting a unique filter word and then tweeting that word from your Twitter account. You will see your ...


-2

Welcome to gis.stackexchange.com! You should have read the guide on asking questions: Have you thoroughly searched for an answer before asking your question? Sharing your research helps everyone. Tell us what you found and why it didn’t meet your needs. This demonstrates that you’ve taken the time to try to help yourself, it saves us from ...


2

The actual problem is to plot a circle with 40 miles in diameter on a map with a lat/lon projection (typically EPSG:4326), because native map units are degrees. Therefore, it seems to me the simplest solution is to work with a different projection that is based on meters (that can easily be converted to miles) rather than degrees. As an alternative to ...


1

Here is a very simple and generic example. You need to access the relevant slots of the individual polygons. If you provide a data sample that better fits your situation I can try and be more specific. mp <- readWKT("MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),((15 5, 40 10, 10 20, 5 10, 15 5)))") # total area sapply(slot(mp, "polygons"), slot, "area") ...



Top 50 recent answers are included