Hot answers tagged raster
6
You have to decide whether you want to enter coordinates in lat/lon, or from a background map with a projected CRS (e.g. in metres). For the first choice, you have to set the target CRS to WGS84 (EPSG:4326), and don't mix up lat and long values.
For the second choice, the target CRS must be set to the same CRS as the project CRS. This can be different from ...
5
I will offer an R solution that is coded in a slightly non-R way to illustrate how it might be approached on other platforms.
The concern in R (as well as some other platforms, especially those that favor a functional programming style) is that constantly updating a large array can be very expensive. Instead, then, this algorithm maintains its own private ...
5
Plotting the estimated slopes, as in the question, is a great thing to do. Rather than filtering by significance, though--or in conjunction with it--why not map out some measure of how well each regression fits the data? For this, the mean squared error of the regression is readily interpreted and meaningful.
As an example, the R code below generates a ...
5
I would recommend clipping the raster to the shapefile, then in the resulting raster you can look at the number of cells present for each of your classes. The area can be calculated by multiplying the number of cells by the area covered each pixel (cell size squared).
It's a different approach than the equally valid solution offered above but from a ...
4
Depending on the situation you have a few options:
1) The Region Group tool will allow you to isolate individual pixels, which can then be filtered based on count.
2) The Con tool can be used to remove (set null) or change the value of pixels if you can identify them by their values.
3) A variety of filters can be used to remove isolated pixels (e.g. ...
4
I'd recommend using the GDAL Pytho API as a substitute for ArcPy, which will give you an extensive range of functions. If you are using GDAL then you have quite an extensive list of possible formats but GeoTiff is always a popular choice but there are plenty of other options (check the GDAL documentation).
Alternatively, if you don't want to use GDAL, then ...
4
You could use a masked array like this:
#Assuming your NoData value is −32768
aa=numpy.ma.MaskedArray(arcpy.RasterToNumPyArray("water_depth"))
aa.mask=(aa==-32768)
#OR
#Explicitly set the nodata value
aa=numpy.ma.MaskedArray(arcpy.RasterToNumPyArray("water_depth",nodata_to_value=-999))
aa.mask=(aa==-999)
#OR
#Or just mask anything negative
...
3
You can carry this out quickly using the Raster Calculator, either using the ArcMap user interface or a python script (see previous link for the code snippet).
You already know your cell resolution (30m by 30m) and your depth (represented by each cell in the raster), so you simply need to calculate the volume of each pixel as a rectangular prism. Make sure ...
3
In 2003 the Landsat 7 ETM sensor had an error in the scan line corrector. The gaps you are observing are a result of this sensor error. There are many ways to correct this issue. Essentially it requires that you find two similar images and fill in the gaps. In other words, there is another preprocessing step that is required. Some resources that may ...
3
This is a concise way to do that in R --- here without intermediate files:
library(raster)
raster_data <- list.files(path=getwd()) #promt user for dir containing raster files
s <- stack(raster_data)
f <- function(x) { rowSums(x >= 4 & x <= 9) }
x <- calc(s, f, progress='text', filename='output.tif')
3
Assuming that presencias and variables share the same projection, this should be an easy task. I recommend you to add these lines of code after your read.table() statement in order to convert presencias dataframe to a SpatialPointsDataFrame object (just refine the names of the columns containing x and y coordinates if they differ from my example).
...
3
Insert 'import arcpy' at the top of your code (NameError: name 'arcpy' is not defined)
No indentation needed after the import statements
Check out the ArcGIS Spatial Analyst extension license using 'arcpy.CheckOutExtension("Spatial")' (RuntimeError: ERROR 000824: The tool is not licensed.)
Assuming that ...
3
This can be accomplished using three main steps:
Convert your raster to polygons using the Raster To Polygon tool
Overlay your polygons (from raster) with your shapefile using the Intersect tool or Union tool - which depends on how many polygons you want output
Tabulate your results using the Summary Statistics tool (or perhaps the Tabulate Intersections ...
3
If your raster has cells of equal size then the area is easy enough, as its just the cells size squared. Then you could use the Spatial Analyst - Math - Times tool to multiply each cell times its area to get a raster representing the volume for each cell.
import arcpy
from arcpy.sa import *
outVol = Times("depth", "100")
One method to get the total would ...
3
A profile of the horizon plots the apparent elevation of the land-sky demarcation against the direction of view (the "azimuth").
In this plot the "adjusted altitude" measures the angle of view (shown as 1000 times its tangent). It was obtained from a DEM by first computing the viewshed for a 20 meter fire tower at a location near the middle:
This ...
3
You can encode multiple values in raster data using bit flags. This is how MODIS data stores quality and cloud masks. Bitwise operators can be used to extract these values.
3
I suspect much of what you are seeing as inconsistent software behaviour will become clearer after reviewing the online help for Raster dataset attribute tables. The help for Build Raster Attribute Table is probably also worth reviewing.
If the solution does not become evident, then I would recommend listing the precise steps you performed on a raster ...
3
I would drop using the maps package and find a state shapefile. Then load that into R using rgdal, and then do some polygon overlay work.
# use state bounds from gadm website:
us=shapefile("USA_adm1.shp")
# extract states (need to uppercase everything)
ne = us[match(toupper(nestates),toupper(us$NAME_1)),]
# create a random raster over the space:
...
2
I ended up using the ArcGIS floating-point raster format (extension .flt), as described here.
The C++ code I wrote to handle things follows. It will not compile without a grid class, but it should be easy enough to read and understand.
/**
@brief Writes a floating-point grid file
@param[in] &basename Name, without extension, of output file
...
2
Adding my comment as an answer since this seems to be a solution:
ArcMap 10 sometimes has issues with background processing on certain machines. If you have background processing enabled, try disabling it (Geoprocessing -> Geoprocessing Options -> Uncheck Enable under Background Processing). Try running your process again and see if that works.
2
An alternate approach to the proposed DEM extraction, assuming you have Spatial Analyst, would be to:
Use the 'Slice' tool with the "EQUAL AREA" method and 5 or 10 zones.
Then you could either:
a) Re-classify the Slice output to create a mask to extract your DEM
-- or --
b) Use it with the 'CON' tool or Raster Calculator to mask out your DEM on the top ...
2
Using Spatial Analyst, we can expect each iteration of IDW, thresholding, and accumulation of results to take from a fraction of a second to a few seconds. Multiplying by 24 * 365 = 8760 hours in a year suggests this calculation will take hours at the least.
It can be much faster. Faster is not always better, but here--because there are significant issues ...
2
I'm not sure about Tile Mill, but if you have ArcMap you can do this fairly easily by importing the image and Geo-referencing it, the advantage of this is if you're using version 10 or later you can use ArcMap's base data feature to use satellite imagery as a background, and then just align your floor plan based on that.
Here's a guide on how to ...
2
As I noticed in comments, generally you should avoid using R for non-statistical purposes due to performance issues in certain aspects (working with cycles is an example). Here is code example for you in Pyhton (thanks to this article) for reclassification of a single file with a single band. You will be able to modify it easily for batch processing if you ...
2
DEMs have 1 band which hold elevation values. For comparison, a color image has three bands (Red, Green, Blue). You can check information about any raster dataset by right-clicking on the layer > Layer Properties > Source. Other particularly useful information located here include pixel depth, cell size, coordinate system, and format. Make sure to set ...
2
Have you tried inspecting the page's source (e.g. with Firefox's developer tools)? Looks like it's using OpenLayers and GeoExt. The time slider appears to be a custom development (check out the functions in the GeoAdmin namespace in the app.js source -- you will probably want to use something like jsbeautifier to make it readable).
2
You can do
library(raster)
r <- raster(matrix)
# replace with correct coordinates
extent(r) <- (0, 1, 0, 1)
r <- writeRaster(r, 'filename.tif')
But your question suggests that you would have been better of accessing the ncdf file as a RasterLayer object (and avoid creating matrices and keeping track of coordinates)
pft <- ...
2
Sadly, no. Not in 2.0. 2.1 should provide significant performance improvement over 2.0.
No matter what you attempt in the database with 2.0, you are limited by the number of times the raster objects are serialized/deserialized.
PostGIS 2.1 should be going beta in the next week or two (hopefully sooner than later).
2
Without ever having performed your operation, and no time to spare for play, I can only add these two links to your list:
Find the Nearest Raster Cell Value Based on a Vector Point (The first answer (with 4 votes) is what intrigued me).
Also: would Hawth's Gridspread help?
2
By making a buffer you approximate the distances (as you do not get a true circle, but a polygonal approximation of one). You can improve the rgeos answer by increasing the number of segments of the circle with the "quadsegs" argument:
nbg.bff2 <- gBuffer(nuremberg, quadsegs=50, width = 20000)
table( unlist(extract(germany, nbg.bff2)) )
0 1
118 91 ...
Only top voted, non community-wiki answers of a minimum length are eligible
