Here's a curly one!  EDIT: As an alternative the modal raster value inside each bounding box would work, and might possibly be faster. Still unsure how to do it though. Perhaps Zonal statistics using Majority. Then I'll need to filter by the maximum count.   Background Info Currently using ArcGIS Pro.     I currently have a grid, and a map of vegetation in raster format. I want to find the most common vegetation by area, in each grid square. The dataset is probably a couple of hundred thousand rows long - and likely isn't nicely sorted.    Current Process Vegetation Raster to Vegetation Polygons (Raster to Polygon) Pairwise Intersect Vegetation Polygons with Grid Polygons to create VegetationGrid_Polygons Calculate Geometry Attributes on VegetationGrid_Polygons to give Area Stuck.  Starting Table Now I should theoretically have a table that looks like: VegetationID, GridID, Area 1, 10, 2 2, 10, 8 &lt;-- Notice this, there are two VegetationID's #2 2, 10, 4 &lt;-- 5, 10, 4 1, 5, 10 3, 5, 46 5, 5, 4    What tool(s) can I use to create a table that has the sum area of each VegetationID by each GridID? VegetationID, GridID, Area 1, 10, 2 2, 10, 12 &lt;-- This is now 12 5, 10, 4 1, 5, 10 3, 5, 46 5, 5, 4    Then, how to find the VegetationID with the largest area value by GridID? VegetationID, GridID, Area 2, 10, 12 3, 5, 46
15101
First, I have a shapefile consisting of square grids and each grid has an ID. Second, I have a shapefile consisting of different zones (with IDs).  A grid ID in the first shapefile can contain many zone IDs from the second shapefile.  I want to know what is the dominant zone ID inside a grid and record that in a table. I only want to know the dominant zone (in terms of the area) and not the rest of the zones intersecting that grid. In the end, I would want a 2-column table with all the grid IDs in the first column and the corresponding dominant zone ID (or the largest zone in terms of area) in the second column.  Here is a sample figure.  By just looking at the figure, I could tell that in Grid ID 1, the dominant zone/area is Zone ID 1.  In Grid ID 2, the dominant zone/area is Zone ID 2.  And in Grid ID 3, the dominant zone/area is Zone ID 3.  I would like to have a table with two columns, with each grid showing only the ID of the dominant area.
15934
I've a set of regions within a city. I'd like to split each region into five subregions and then get the centroid of each one of them.  I know how to make a vector grid, but the problem I have is that the generated points are not taking into account the regions, they're only related to the general map.  I have this:    And I want to split each region in five equal size areas. Something like this:    I'd like to build a visualization similar to this one: . In order to accomplish that, I need to split each region, because the original ones are pretty big.
97327
I am using ArcGIS 9.3 and python 2.5. My code is shown at the bottom of this post.  I have a raster containing fire intensity information at 200m cell resolution. I run the script and the first count cell part returns the size of the fire. This works fine.  I then clip the raster using a polygon representing a forested catchment. The new raster represents the area of the catchment that is burnt. When I try to use the same code to count the number of cells of the new raster, it always gives me an answer of zero, even if there is definitely a raster created by the clip tool.  I have tried every combination of every tool I can think of.  my code...      fire = str(ignitionID) + "_" + str(weatherID) + "_Intensity.asc"     print "fire: " + fire      # ---------------------------------------------------------------------------     # this section calculates the size of the fire for model calibration purposes     # ---------------------------------------------------------------------------      # Process: Build Raster Attribute Table...     the_fire = input_fires + fire     shutil.copy(the_fire, "C:\\FIRESTORM\\Workspace") # copy the fire to workspace     the_fire = "C:\\FIRESTORM\\Workspace\\" + fire      try:       gp.BuildRasterAttributeTable_management(the_fire, "NONE")       resultf = gp.GetCount_management(the_fire)  # count num of rows in raster       countf = int(resultf.GetOutput(0))       fire_area = countf * 4  # each cell has an area of 4 ha       print "the fire area = " + str(fire_area) + " ha"      except:  # the attribute table cannot be created coz the fire is too large        fire_area = 300000 # default value for large fires 300000 ha       print "the fire area EXCEEDS 262,140 ha"      # -----------------------------------------------------------------------------     # this section calculates the area of the catchments that are burnt by the fire     # -----------------------------------------------------------------------------      # Local variables...     catch_extent = input_catch + "WS_Catchment_region.shp"  # catchment boundary     fire_raster = "C:\\FIRESTORM\\Workspace\\fire_raster"     fire_ascii = the_fire     Output_raster = "C:\\FIRESTORM\\Workspace\\Extract_ASCI1"     Output_ASCII = "C:\\FIRESTORM\\Workspace\\rastert_extract1.asc"     catch_burn_area = 0  # reset catchment burn area      # Process: ASCII to Raster...     gp.ASCIIToRaster_conversion(fire_ascii, fire_raster, "INTEGER")      # Process: Define Projection... based on projection of catch_extent     desc = gp.Describe(catch_extent)     spatialRef = desc.SpatialReference     gp.DefineProjection_management(fire_raster, spatialRef)      # Process: Extract by Mask...     gp.ExtractByMask_sa(fire_raster, catch_extent, Output_raster)      # Convert raster to ascii     gp.RasterToASCII_conversion(Output_raster, Output_ASCII)      # Process: Build Raster Attribute Table...     gp.BuildRasterAttributeTable_management(Output_ASCII, "OVERWRITE")      # count number of rows in raster     resulti = gp.GetCount_management(Output_ASCII)     counti = int(resulti.GetOutput(0))     catch_burn_area = counti * 4  # each cell has an area of 4 ha      print "catchment area burnt = " + str(catch_burn_area) + " ha"
110044
I have downloaded raster map data from ONS, the data has come in the form of two files... one is full of  .tlw files and the other is full of .tab files...   How do I import it into a raster layer to ArcCatalog?
102902
I have two layers. A polygon-shape-layer with many tiles and a raster-layer containing  with many categories in a colourmap. I want to obtain for each polygon in the shapelayer a sum of each land cover category of the raster-layer.  For example there is a polygon with id '2' and i want to Attributes like this for this polygon (in percent or square meters):   Arable land : 15 % Forest: 11 % Streets:2 % (... and so one)   I tried to do it in grass, qgis (no function), saga (just sums up every to a total value) r(total sum), but i still found no solution. Most plugins (zonal statistics in qgis) only support 0-1 raster layers. v.rast.stats didn't help either. Iam open to any good and smart solution!. Maybe i even used a wrong approach or made mistakes.   In Arcgis this task is quite easy, if am remember right, but i am still missing a good solution for your everyday linux user.  I am running a debian linux system and this why i can only use programs for this OS.    EDIT: Because this question still has so many views and visitors: I wrote a QGIS-plugin, which also is capable of calculating the landcover of raster layer. I have'nt coded a polygon overlay yet, but it definitely planed. Find the plugin  and install the Scipy library first.
13255
I have two polygons. One contains fields(X,Y,Z) and the other contains soil types (A,B,C,D). I want to know what  area of every field contains which type of soil. I tried the following:    library(rgdal) library(rgeos) Field&lt;-readOGR("./","Field") Soil&lt;-readOGR("./","Soil") Results&lt;-gIntersects(Soil,Field,byid=TRUE) rownames(Results)&lt;-Field@data$FieldName colnames(Results)&lt;-Soil@data$SoilType  &gt; Results       A     B     C     D Z  TRUE FALSE FALSE FALSE Y FALSE  TRUE  TRUE FALSE X  TRUE  TRUE  TRUE  TRUE   and achieved good results with it telling me which field contains which soil type. However, how do I get the area instead?
110216
In R, it's relatively trivial to perform per pixel calculations based on a raster stack (e.g., get std.dev for each pixel on a 12 layer GeoTIFF). Unfortunately, the speed is less than desirable when working with large rasters.  Recently, I've been trying to use GDAL for this purpose, but I can't seem to get the expression down. I've tried the following, but it doesn't seem to work as advertised (or how I think it should work):  gdal_calc.py -A tmean_monthly_1980.tif --calc="std(A)" --outfile=tstdev_1980.tif --allBands=A --calc="std(A)" --NoDataValue=-9999   In the above example tmean_monthly_1980.tif is a raster composed of 12 layers (Jan-Dec) of mean monthly temperatures in 1980 (stacked using gdal_merge). Although the GDAL command produces an output, the results are definitely not the standard deviation per pixel.  Another option might be to use single band files, but gdal_calc.py appears limited to "only" 26 inputs (A-Z; i.e., "-A file -B file -C file etc."), and the resulting command line/script would be fairly cumbersome when working with large, multi-layered rasters (e.g., time series).  Is there a specific syntax for performing per-pixel raster calculations on a multi-band image using GDAL? If only individual files/bands can be specified using gdal_calc.py, please provide a scripted example using Bash, Python, etc.  Specifically, my goal is to calculate the standard deviation and cv for each pixel in a multi-layer raster stack, or using individual files that make up the same raster stack.  Note: I'm holding out for a GDAL specific answer, but other FOSS solutions are welcomed and will be upvoted, especially if they are relatively compact and adaptable to large number of inputs.
46006
In my daily work, I am constantly asked to calculate areas of global raster datasets in geographic projection at 30 arc second resolution. These datasets are normally the result of a Combine operation (a typical example is a vegetation classes combined with country layer). To do this, our unit created a raster dataset with the area of each pixel in geographic projection at 30 arc second. With this area grid, a zonalstat is performed to sum the areas for each class. Since I am not sure how this area grid was created, I always wondered whether this approach is more accurate of just reprojecting the raster in a equal-area projection (from simple tests the results of the two methods are similar). Does anybody experienced a similar situation?
110627
I have a large number of rasters with 3 bands (R,G,B). What I want to do is add all 3 bands of each raster together to end up with a single band image with the summed R,G,B values in each pixel but I'm not sure how to do this in ModelBuilder. Do I have to somehow split each raster first, then add each one?  How would I iterate through each band of every raster?  Edit:Basically this is what I have so far. Right now it will take 3 bands of 1 raster and add the bands then go on to create a polygon of a certain reclassed area. I want to do this on a large number of rasters but am unsure of how to make it work.
90058
How do I find the coordinates of the centre / centroid of a raster?   I have tried Zonal Geometry but this does not give me the coordinates.  I need to find the centroids of two rasters and create points and measure the distance and bearings between them. i.e. distance and direction between two raster centroids.  I am using ArcGIS Desktop 10.2.
14161
I plan on rasterizing (terminology?) census block data relating to density to a regular 30 m by 30 m grid for use with NLCD data.  I then hope to distribute the density of the original census block among each result set of grid cells within the original polygon according to the land cover type dictated by the NLCD dataset.  For example, let's imagine a string of city blocks with a park, all found within one census block.  If within this census block the average density is 1000 people/km^2, we ultimately would want the resulting grid cell containing the park to have a lower population density than the grid cells containing the surrounding apartments.  This information may be captured by the NLCD dataset, with the park being classified as '21' (open but developed land) versus 24 for apartment buildings which could be classified as '24' (high density development).  Is there a method, or set of methods, within ArcMap that would allow me to:   Redistribute the density based on the NLCD value Perform this redistribution only within the original census block polygon   I don't know if this is possible, but I figured I would ask.  Edit  Here's what I have been trying:  I converted the census block data to a raster, snapping it to the NLCD land cover data and ensuring the output grid cells were the same size.  Each cell contains the same info as the original census block (housing umbers, population, etc). I have been trying to figure out a way to edit the new grid so that if the NLCD grid cell contains a land cover type that doesn't support a business or housing (park, forest, farm, etc.) that it's housing and population values area set to 0.  My next move was going to be figuring out how to subtract the non-residential/non-business grid cells area from the area of the original census block, and force a recalculation of the density.  Converted the altered census block layer to grid cells again, with each having the new density value, and once again set the cells' housing, population, and density values that don't support a business or housing (park, forest, farm, etc.) to 0.  Kind of circular, and I'm not sure of all the steps (raster calculator for editing individual grid cells based on another raster's grid cell value?), but it seems to logically work out if I can get the parts to work.  On the other hand, @dof1985's answer seems to simplify the first part of my workflow (by directing me to the correct toolset), but I will still have to nullify the grid cells where parks, etc., do exits after the density has been recalculated for the census data.
58334
I have a single band raster file as drawn in the labeled image attached.  The size of raster is around 4k pixel in each x and y axis. Large part of this raster is empty. Only in small portions of the raster contains data as labeled. Is there any efficient way I can extract rectangular patches around the data portions of the large tif raster file? Any pythonic way using gdal or rasterio or some other package. What I have done so far: I have followed from sklearn.feature_extraction.imageway to extract 2D patches where there is data as follows. from sklearn.feature_extraction.image import extract_patches_2d im = imageio.imread(largeFilePath) # reading the large image x, y = np.nonzero(im) # taking xy locations with non zero values  # determining xleft, xright, yleft and yright of the non zero area in  raster  xl,xr = x.min(),x.max()  yl,yr = y.min(),y.max()   #Extraing the non zero area from my original raster only_desired_area = im[xl:xr+1, yl:yr+1]  # my desired window to extract  window_shape = (480, 480)  # 2D Ptaches extraction B = extract_patches_2d(only_desired_area, window_shape, max_patches=100)  Issue: This doesn't have georeference information because I haven't used gdal or rasterio
81944
I'm trying to reclassify only a certain part of a raster...not the entire raster.  Here's the problem.  I have a LULC raster of the entire US.  There are some issues with the raster that I want to fix, namely some cells that should be water are classified as wetlands.  These are mainly off the coast, not inland.  However, some cells that are classified as wetlands should really be wetlands (mainly inland) and I don't want these to change.  I plan on using a ZIP Code boundary file to determine which cells are considered "off the coast" and which ones are not.  I have access to ArcGIS 10.1 and R with the appropriate raster packages I need so a solution in either software would be fine.  Can anyone help me figure this out?
9278
I have a raster (USGS DEM actually) and I need to split it up into smaller chunks like the image below shows. That was accomplished in ArcGIS 10.0 using the Split Raster tool. I would like a FOSS method to do this. I've looked at GDAL, thinking surely it would do it (somehow with gdal_translate), but can't find anything. Ultimately, I'd like to be able to take the raster and say how large (4KM by 4KM chunks) I would like it split up into.
78049
I found  in the old ESRI forums with the same type of problem I am facing today. The solution to the problem is explained in little detail and since I am not an expert, I was wondering if somebody could describe it to me in little more detail.  We are doing a mental map analysis. We have let's say 50 polygons of a certain city square which are all different in shape, size and location. Our goal is to create a single polygon which would represent the answers from our whole input data. So basically, we need a method for generating a single polygon out of these 50 which would be an "average" polygon in size, shape and location.  Can somebody explain the process?
33271
I am looking for a method to process a remote sensing image and extract the crown areas of the individual trees from the image.   I have both visual wavelength areal imagery, and lidar data from the area. The location in question is a desert area, so the tree cover isn't as dense as a forest area. The resolution of the aerial imagery is 0.5 feet by 0.5 feet. The lidar resolution is approximately 1 x 1 feet. Both the visual data and the lidar come from a Pima County, Arizona dataset. A sample of the type of aerial imagery I have is at the end of this post.   This question  seems to be the same issue, but there does not seem to be a good answer there.  I can obtain a reasonable classification of the vegetation types (and information about the overall percent cover) in the area by using the Iso Cluster classification in Arcmap, but this provides little information on individual trees. The closest I have to what I want is the results of passing the output of the isocluster classification through the Raster to Polygon feature in Arcmap. The problem is that this method merges near by trees into a single polygon.   Edit: I probably should have included some more detail about what I have. The raw datasets I have are:    Full las data, and a tiff raster generated from it. Visual imagery (like the sample image shown, but covering a much wider area)  Manual direct measurements of a subset of the trees in the area.    From these I have generated:    The ground/vegetation classifications. The DEM/DSM rasters.
19042
I have a large collection of points with X and Y coordinates in a projected coordinate reference system (EPSG:28992). I would like to compute a square buffer around it (say, 500 meters around the point), so I can crop this polygon with a raster layer (in TIFF format) and proceed with a further analysis with the values within the cropped windows.  I am aware that this is possible in , but for the sake of the workflow I am implementing, I prefer to have it entirely in the Python GDAL environment, so it becomes automatic. Since I know the X/Y coordinates and the size of the grid cells in the raster, I could add/subtract to the coordinates the desired buffer, but I wonder whether there is a less handcrafted version of this "square-buffer" operation, which is less prone to errors. Other Python libraries (e.g. shapely) would be fine for me as well.   Note that a band.GetRasterBand(1).ReadAsArray(X, Y, 5, 5) does not return the window around the point, but the window leaving the point in a corner.   For instance, the following code:  path_curr_lyr = r"path_to_tif_file" tif = gdal.Open(path_curr_lyr) dat = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 1, 1) print(dat) sq_buf = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 5, 5) print(sq_buf)   Accesses the raster grid cell X/Y (9550, 9620) and returning the value of forest coverage in that pixel, and the value of a window of width= 5 pixels.  &gt;&gt;&gt; [[44]] &gt;&gt;&gt; [[ 44  87  99  86  43]      [ 96  99  99  79  14]      [ 99  99  98  34   9]      [ 99  86  27  93  71]      [ 98  98 100  97  65]]    As you can see, this approach returns a window in which the "44" is not a central value, but one in the corner. I would like to have a procedure in which a window around the given grid cell coordinate is created.  Any advice on how to implement this?