Tag Info

Hot answers tagged

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 ...


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 ...


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 ...


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 ...


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 ...


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"))


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

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 ...


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 ...


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

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.


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 ...


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

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 ...


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 ...


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 <- ...


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.


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 ...


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

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. ...


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 ...


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 <- ...


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 ...


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

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

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/


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 <- ...


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 ...


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 ...


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") ...



Only top voted, non community-wiki answers of a minimum length are eligible