Finding minimum bounding extent of given pixel value within raster? The Next CEO of Stack OverflowGetting a raster extent excluding nodata areaAutomatically discard nodata area of rasterCreating Latitude Grid from DEMIs there a “trim” function in Gdal which would remove a NODATA borderRaster Attribute Table Where Each Pixel = 1 RowCropping or removing borders (not NoData) from multiple scanned and georeferenced maps (ArcMap or QGIS)Assign nearest raster cell value to pointHow do I change the extent of a raster file in ArcGIS in order to remove nodata values?Converting a raster to shapefile for only one value?Extracting raster by mask outputs raster with extent beyond the mask layer's in QGISCalculating area of raster pixel when pixels are different sizes at different latitudes using ArcGIS Desktop?Cell value changes during raster clipping using QGIS?Reclassifying raster/extracting pixel values conditionall over multipolygon shapefileExtracting a single value to a new raster from a singleband pseudocolor raster image

What does "Its cash flow is deeply negative" mean?

Has this building technique been used in an official set?

How to count occurrences of text in a file?

Complex fractions

sp_blitzCache results Memory grants

Inappropriate reference requests from Journal reviewers

Is it my responsibility to learn a new technology in my own time my employer wants to implement?

How to start emacs in "nothing" mode (`fundamental-mode`)

Do I need to enable Dev Hub in my PROD Org?

Should I tutor a student who I know has cheated on their homework?

Phase of a real number

Why is the US ranked as #45 in Press Freedom ratings, despite its extremely permissive free speech laws?

In excess I'm lethal

Why do professional authors make "consistency" mistakes? And how to avoid them?

Would a completely good Muggle be able to use a wand?

Why am I allowed to create multiple unique pointers from a single object?

Why don't programming languages automatically manage the synchronous/asynchronous problem?

Is micro rebar a better way to reinforce concrete than rebar?

MessageLevel in QGIS3

How to solve a differential equation with a term to a power?

Sending manuscript to multiple publishers

Why has the US not been more assertive in confronting Russia in recent years?

If a black hole is created from light, can this black hole then move at speed of light?

Does it take more energy to get to Venus or to Mars?



Finding minimum bounding extent of given pixel value within raster?



The Next CEO of Stack OverflowGetting a raster extent excluding nodata areaAutomatically discard nodata area of rasterCreating Latitude Grid from DEMIs there a “trim” function in Gdal which would remove a NODATA borderRaster Attribute Table Where Each Pixel = 1 RowCropping or removing borders (not NoData) from multiple scanned and georeferenced maps (ArcMap or QGIS)Assign nearest raster cell value to pointHow do I change the extent of a raster file in ArcGIS in order to remove nodata values?Converting a raster to shapefile for only one value?Extracting raster by mask outputs raster with extent beyond the mask layer's in QGISCalculating area of raster pixel when pixels are different sizes at different latitudes using ArcGIS Desktop?Cell value changes during raster clipping using QGIS?Reclassifying raster/extracting pixel values conditionall over multipolygon shapefileExtracting a single value to a new raster from a singleband pseudocolor raster image










9















I am wondering if there is a way to find the minimum bounding extent for a raster with a particular value. I clipped a raster from a global image and the extent is set as the global extent with a lot of NoData area. I would like to remove the NoData area from this raster and retain only the area containing the pixels of the particular value. How can I do this?



Here is my example: I would like to extract value = 1 (Blue area) and use the blue area's extent rather than the whole world for further processing.



Sample Image










share|improve this question
























  • Could you post a sample?

    – Aaron
    Jan 4 '13 at 0:36











  • "I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

    – blah238
    Jan 4 '13 at 0:54











  • By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

    – blah238
    Jan 4 '13 at 0:55







  • 1





    @Tomek, the OP is looking to find the extent, not have to create it manually.

    – blah238
    Jan 4 '13 at 9:40







  • 1





    If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

    – whuber
    Mar 21 '13 at 21:00















9















I am wondering if there is a way to find the minimum bounding extent for a raster with a particular value. I clipped a raster from a global image and the extent is set as the global extent with a lot of NoData area. I would like to remove the NoData area from this raster and retain only the area containing the pixels of the particular value. How can I do this?



Here is my example: I would like to extract value = 1 (Blue area) and use the blue area's extent rather than the whole world for further processing.



Sample Image










share|improve this question
























  • Could you post a sample?

    – Aaron
    Jan 4 '13 at 0:36











  • "I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

    – blah238
    Jan 4 '13 at 0:54











  • By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

    – blah238
    Jan 4 '13 at 0:55







  • 1





    @Tomek, the OP is looking to find the extent, not have to create it manually.

    – blah238
    Jan 4 '13 at 9:40







  • 1





    If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

    – whuber
    Mar 21 '13 at 21:00













9












9








9


4






I am wondering if there is a way to find the minimum bounding extent for a raster with a particular value. I clipped a raster from a global image and the extent is set as the global extent with a lot of NoData area. I would like to remove the NoData area from this raster and retain only the area containing the pixels of the particular value. How can I do this?



Here is my example: I would like to extract value = 1 (Blue area) and use the blue area's extent rather than the whole world for further processing.



Sample Image










share|improve this question
















I am wondering if there is a way to find the minimum bounding extent for a raster with a particular value. I clipped a raster from a global image and the extent is set as the global extent with a lot of NoData area. I would like to remove the NoData area from this raster and retain only the area containing the pixels of the particular value. How can I do this?



Here is my example: I would like to extract value = 1 (Blue area) and use the blue area's extent rather than the whole world for further processing.



Sample Image







arcgis-desktop raster extents






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited Jun 30 '14 at 3:46









PolyGeo

53.8k1781245




53.8k1781245










asked Jan 3 '13 at 23:11









SeenSeen

1,13031833




1,13031833












  • Could you post a sample?

    – Aaron
    Jan 4 '13 at 0:36











  • "I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

    – blah238
    Jan 4 '13 at 0:54











  • By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

    – blah238
    Jan 4 '13 at 0:55







  • 1





    @Tomek, the OP is looking to find the extent, not have to create it manually.

    – blah238
    Jan 4 '13 at 9:40







  • 1





    If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

    – whuber
    Mar 21 '13 at 21:00

















  • Could you post a sample?

    – Aaron
    Jan 4 '13 at 0:36











  • "I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

    – blah238
    Jan 4 '13 at 0:54











  • By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

    – blah238
    Jan 4 '13 at 0:55







  • 1





    @Tomek, the OP is looking to find the extent, not have to create it manually.

    – blah238
    Jan 4 '13 at 9:40







  • 1





    If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

    – whuber
    Mar 21 '13 at 21:00
















Could you post a sample?

– Aaron
Jan 4 '13 at 0:36





Could you post a sample?

– Aaron
Jan 4 '13 at 0:36













"I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

– blah238
Jan 4 '13 at 0:54





"I would like to delete the null rows and columns for this raster." What exactly does this mean? I don't understand what the desired end result is.

– blah238
Jan 4 '13 at 0:54













By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

– blah238
Jan 4 '13 at 0:55






By "minimum bounding extent" are you looking for a rectangular extent or a polygonal footprint representing the area of the image with data?

– blah238
Jan 4 '13 at 0:55





1




1





@Tomek, the OP is looking to find the extent, not have to create it manually.

– blah238
Jan 4 '13 at 9:40






@Tomek, the OP is looking to find the extent, not have to create it manually.

– blah238
Jan 4 '13 at 9:40





1




1





If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

– whuber
Mar 21 '13 at 21:00





If literally anything is fair game, then some software has built-in commands to do this; see reference.wolfram.com/mathematica/ref/ImageCrop.html for instance.

– whuber
Mar 21 '13 at 21:00










6 Answers
6






active

oldest

votes


















6














IF I have understood the question correctly it sounds like you want know the minimum bounding box of the values that are not null.Maybe you could convert the raster to polygons, select the polygons you are interested in and then convert them back to a raster. You can then look at the properties values which should give you the minium bounding box.






share|improve this answer


















  • 1





    All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

    – blah238
    Jan 4 '13 at 9:35











  • I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

    – Seen
    Jan 4 '13 at 15:25











  • are you after a python solution?

    – dango
    Jan 7 '13 at 11:52



















8














The trick is to compute the limits of the data that have values. Perhaps the fastest, most natural, and most general way to obtain these is with zonal summaries: by using all non-NoData cells for the zone, the zonal min and max of grids containing the X and Y coordinates will provide the full extent.



ESRI keeps changing the ways in which these calculations can be done; for instance, built-in expressions for the coordinate grids were dropped with ArcGIS 8 and appear not to have returned. Just for fun, here's a set of fast, simple calculations that will do the job no matter what.




  1. Convert the grid into a single zone by equating it with itself, as in



    "My grid" == "My grid"



  2. Create a column index grid by flow-accumulating a constant grid with value 1. (The indexes will start with 0.) If desired, multiply this by the cellsize and add the x-coordinate of the origin to obtain an x-coordinate grid "X" (shown below).


  3. Similarly, create a row index grid (and then a y-coordinate grid "Y") by flow-accumulating a constant grid with value 64.


  4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y": you now have the desired extent.


Final image



(The extent, as shown in the two tables of zonal statistics, is depicted by a rectangular outline in this figure. Grid "I" is the zone grid obtained in step (1).)



To go further, you will need to extract these four numbers from their output tables and use them to limit the analysis extent. Copying the original grid, with the restricted analysis extent in place, completes the task.






share|improve this answer






























    8














    Here's a version of @whubers method for ArcGIS 10.1+ as a python toolbox (.pyt).



    import arcpy

    class Toolbox(object):
    def __init__(self):
    """Define the toolbox (the name of the toolbox is the name of the
    .pyt file)."""
    self.label = "Raster Toolbox"
    self.alias = ""

    # List of tool classes associated with this toolbox
    self.tools = [ClipNoData]


    class ClipNoData(object):
    def __init__(self):
    """Clip raster extent to the data that have values"""
    self.label = "Clip NoData"
    self.description = "Clip raster extent to the data that have values. "
    self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

    self.canRunInBackground = False

    def getParameterInfo(self):
    """Define parameter definitions"""
    params = []

    # First parameter
    params+=[arcpy.Parameter(
    displayName="Input Raster",
    name="in_raster",
    datatype='GPRasterLayer',
    parameterType="Required",
    direction="Input")
    ]

    # Second parameter
    params+=[arcpy.Parameter(
    displayName="Output Raster",
    name="out_raster",
    datatype="DERasterDataset",
    parameterType="Required",
    direction="Output")
    ]

    return params

    def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return arcpy.CheckExtension('spatial')==u'Available'

    def execute(self, parameters, messages):
    """See https://gis.stackexchange.com/a/55150/2856
    ##Code comments paraphrased from @whubers GIS StackExchange answer
    """
    try:
    #Setup
    arcpy.CheckOutExtension('spatial')
    from arcpy.sa import *
    in_raster = parameters[0].valueAsText
    out_raster = parameters[1].valueAsText

    dsc=arcpy.Describe(in_raster)
    xmin=dsc.extent.XMin
    ymin=dsc.extent.YMin
    mx=dsc.meanCellWidth
    my=dsc.meanCellHeight
    arcpy.env.extent=in_raster
    arcpy.env.cellSize=in_raster
    arcpy.AddMessage(out_raster)

    ## 1. Convert the grid into a single zone by equating it with itself
    arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
    zones = Raster(in_raster) == Raster(in_raster)

    ## 2. Create a column index grid by flow-accumulating a constant grid
    ## with value 1. (The indexes will start with 0.) Multiply this by
    ## the cellsize and add the x-coordinate of the origin to obtain
    ## an x-coordinate grid.
    arcpy.AddMessage(r'Create a constant grid...')
    const = CreateConstantRaster(1)

    arcpy.AddMessage(r'2. Create an x-coordinate grid...')
    xmap = (FlowAccumulation(const)) * mx + xmin

    ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
    ## constant grid with value 64.
    arcpy.AddMessage(r'3. Create a y-coordinate grid...')
    ymap = (FlowAccumulation(const * 64)) * my + ymin

    ## 4. Use the zone grid from step (1) to compute the zonal min and
    ## max of "X" and "Y"
    arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

    xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORYxrange", "NODATA", "MIN_MAX")
    xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORYxrange", ["MIN","MAX"]).next()

    yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORYyrange", "NODATA", "MIN_MAX")
    ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORYyrange", ["MIN","MAX"]).next()

    arcpy.Delete_management(r"IN_MEMORYxrange")
    arcpy.Delete_management(r"IN_MEMORYyrange")

    # Write out the reduced raster
    arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
    output = Raster(in_raster) * 1
    output.save(out_raster)

    except:raise
    finally:arcpy.CheckInExtension('spatial')





    share|improve this answer

























    • Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

      – matt wilkie
      Mar 22 '13 at 15:52






    • 1





      I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

      – user2856
      Mar 22 '13 at 20:31











    • good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

      – matt wilkie
      Mar 22 '13 at 20:42











    • +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

      – whuber
      Mar 23 '13 at 3:34


















    4














    I've devised a gdal and numpy based solution. It breaks the raster matrix into rows and columns and drop any empty row/column. In this implementation "empty" is anything less than 1, and only single band rasters accounted for.



    (I realise as I write that this scanline approach is only suitable for images with nodata "collars". If your data are islands in seas of nulls, the space between islands will get dropped as well, squishing everything together and totally messing up the georeferencing.)



    The business parts (needs fleshing out, won't work as is):





     #read raster into a numpy array
    data = np.array(gdal.Open(src_raster).ReadAsArray())
    #scan for data
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
    # assumes data is any value greater than zero
    crop_box = (min(non_empty_rows), max(non_empty_rows),
    min(non_empty_columns), max(non_empty_columns))

    # retrieve source geo reference info
    georef = raster.GetGeoTransform()
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    # write to disk
    band = out_raster.GetRasterBand(1)
    band.WriteArray(new_data)
    band.FlushCache()
    out_raster = None


    In a full script:



    import os
    import sys
    import numpy as np
    from osgeo import gdal

    if len(sys.argv) < 2:
    print 'n [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
    sys.exit(1)

    src_raster = sys.argv[1]
    out_raster = sys.argv[2]

    def main(src_raster):
    raster = gdal.Open(src_raster)

    # Read georeferencing, oriented from top-left
    # ref:GDAL Tutorial, Getting Dataset Information
    georef = raster.GetGeoTransform()
    print 'nSource raster (geo units):'
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]
    cols, rows = raster.RasterYSize, raster.RasterXSize
    print ' Origin (top left): :10, :10'.format(xmin, ymax)
    print ' Pixel size (x,-y): :10, :10'.format(xcell, ycell)
    print ' Columns, rows : :10, :10'.format(cols, rows)

    # Transfer to numpy and scan for data
    # oriented from bottom-left
    data = np.array(raster.ReadAsArray())
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
    crop_box = (min(non_empty_rows), max(non_empty_rows),
    min(non_empty_columns), max(non_empty_columns))

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    new_rows, new_cols = new_data.shape # note: inverted relative to geo units
    #print cropped_transform

    print 'nCrop box (pixel units):', crop_box
    print ' Stripped columns : :10'.format(cols - new_cols)
    print ' Stripped rows : :10'.format(rows - new_rows)

    print 'nCropped raster (geo units):'
    print ' Origin (top left): :10, :10'.format(new_xmin, new_ymax)
    print ' Columns, rows : :10, :10'.format(new_cols, new_rows)

    raster = None
    return new_data, cropped_transform


    def write_raster(template, array, transform, filename):
    '''Create a new raster from an array.

    template = raster dataset to copy projection info from
    array = numpy array of a raster
    transform = geo referencing (x,y origin and pixel dimensions)
    filename = path to output image (will be overwritten)
    '''
    template = gdal.Open(template)
    driver = template.GetDriver()
    rows,cols = array.shape
    out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
    out_raster.SetGeoTransform(transform)
    out_raster.SetProjection(template.GetProjection())
    band = out_raster.GetRasterBand(1)
    band.WriteArray(array)
    band.FlushCache()
    out_raster = None
    template = None

    if __name__ == '__main__':
    cropped_raster, cropped_transform = main(src_raster)
    write_raster(src_raster, cropped_raster, cropped_transform, out_raster)


    The script is in my code stash on Github, if the link goes 404 hunt around a bit; these folders are ripe for some reorganization.






    share|improve this answer




















    • 1





      This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

      – user32882
      Oct 24 '18 at 13:05



















    2














    For all its analytical power, ArcGIS lacks basic raster manipulations that you can find with traditional desktop image editors like GIMP. It expects that you want to use the same analysis extent for your output raster as your input raster, unless you manually override the Output Extent environment setting. Since this is exactly what you are looking to find, not set, the ArcGIS way of doing things is getting in the way.



    Despite my best efforts, and without resorting to programming, I could find no way to get the extent of your desired subset of the image (without raster-to-vector conversion which is computationally wasteful).



    Instead, I used GIMP to select the blue area using the select by color tool and then inverted the selection, hit Delete to remove the rest of the pixels, inverted the selection again, cropped the image to selection, and finally exported it back out to PNG. GIMP saved it as a 1-bit depth image. The result is below:



    Output



    Of course, because your sample image lacked a spatial reference component, and GIMP is not spatially aware, the output image is about as useful as your sample input. You will need to georeference it for it to be of use in a spatial analysis.






    share|improve this answer























    • Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

      – whuber
      Jan 4 '13 at 14:45











    • @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

      – Seen
      Jan 4 '13 at 15:24











    • @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

      – whuber
      Jan 4 '13 at 16:00


















    1














    Here is one possibility using SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021




    In SAGA GIS there is a "Crop to Data" module (in Grid Tools module
    library), which performs the task.



    But this would require you to import your Geotif with the GDAL import
    module, process it in SAGA, and finally export it as Geotif again with
    the GDAL export module.




    Another possibility using only ArcGIS GP tools would be to build a TIN from your raster using Raster to TIN, compute its boundary using TIN Domain, and Clip your raster by the boundary (or its envelope using Feature Envelope to Polygon).






    share|improve this answer

























      Your Answer








      StackExchange.ready(function()
      var channelOptions =
      tags: "".split(" "),
      id: "79"
      ;
      initTagRenderer("".split(" "), "".split(" "), channelOptions);

      StackExchange.using("externalEditor", function()
      // Have to fire editor after snippets, if snippets enabled
      if (StackExchange.settings.snippets.snippetsEnabled)
      StackExchange.using("snippets", function()
      createEditor();
      );

      else
      createEditor();

      );

      function createEditor()
      StackExchange.prepareEditor(
      heartbeatType: 'answer',
      autoActivateHeartbeat: false,
      convertImagesToLinks: false,
      noModals: true,
      showLowRepImageUploadWarning: true,
      reputationToPostImages: null,
      bindNavPrevention: true,
      postfix: "",
      imageUploader:
      brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
      contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
      allowUrls: true
      ,
      onDemand: true,
      discardSelector: ".discard-answer"
      ,immediatelyShowMarkdownHelp:true
      );



      );













      draft saved

      draft discarded


















      StackExchange.ready(
      function ()
      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fgis.stackexchange.com%2fquestions%2f45159%2ffinding-minimum-bounding-extent-of-given-pixel-value-within-raster%23new-answer', 'question_page');

      );

      Post as a guest















      Required, but never shown

























      6 Answers
      6






      active

      oldest

      votes








      6 Answers
      6






      active

      oldest

      votes









      active

      oldest

      votes






      active

      oldest

      votes









      6














      IF I have understood the question correctly it sounds like you want know the minimum bounding box of the values that are not null.Maybe you could convert the raster to polygons, select the polygons you are interested in and then convert them back to a raster. You can then look at the properties values which should give you the minium bounding box.






      share|improve this answer


















      • 1





        All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

        – blah238
        Jan 4 '13 at 9:35











      • I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

        – Seen
        Jan 4 '13 at 15:25











      • are you after a python solution?

        – dango
        Jan 7 '13 at 11:52
















      6














      IF I have understood the question correctly it sounds like you want know the minimum bounding box of the values that are not null.Maybe you could convert the raster to polygons, select the polygons you are interested in and then convert them back to a raster. You can then look at the properties values which should give you the minium bounding box.






      share|improve this answer


















      • 1





        All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

        – blah238
        Jan 4 '13 at 9:35











      • I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

        – Seen
        Jan 4 '13 at 15:25











      • are you after a python solution?

        – dango
        Jan 7 '13 at 11:52














      6












      6








      6







      IF I have understood the question correctly it sounds like you want know the minimum bounding box of the values that are not null.Maybe you could convert the raster to polygons, select the polygons you are interested in and then convert them back to a raster. You can then look at the properties values which should give you the minium bounding box.






      share|improve this answer













      IF I have understood the question correctly it sounds like you want know the minimum bounding box of the values that are not null.Maybe you could convert the raster to polygons, select the polygons you are interested in and then convert them back to a raster. You can then look at the properties values which should give you the minium bounding box.







      share|improve this answer












      share|improve this answer



      share|improve this answer










      answered Jan 4 '13 at 7:12









      dangodango

      654413




      654413







      • 1





        All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

        – blah238
        Jan 4 '13 at 9:35











      • I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

        – Seen
        Jan 4 '13 at 15:25











      • are you after a python solution?

        – dango
        Jan 7 '13 at 11:52













      • 1





        All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

        – blah238
        Jan 4 '13 at 9:35











      • I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

        – Seen
        Jan 4 '13 at 15:25











      • are you after a python solution?

        – dango
        Jan 7 '13 at 11:52








      1




      1





      All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

      – blah238
      Jan 4 '13 at 9:35





      All told this is probably the best approach given the confines of the ArcGIS raster processing workflow.

      – blah238
      Jan 4 '13 at 9:35













      I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

      – Seen
      Jan 4 '13 at 15:25





      I did this exactly. Is there any automatic way? I think the raster to polygon algorithm has a step to extract the minimum bounding box of the raster.

      – Seen
      Jan 4 '13 at 15:25













      are you after a python solution?

      – dango
      Jan 7 '13 at 11:52






      are you after a python solution?

      – dango
      Jan 7 '13 at 11:52














      8














      The trick is to compute the limits of the data that have values. Perhaps the fastest, most natural, and most general way to obtain these is with zonal summaries: by using all non-NoData cells for the zone, the zonal min and max of grids containing the X and Y coordinates will provide the full extent.



      ESRI keeps changing the ways in which these calculations can be done; for instance, built-in expressions for the coordinate grids were dropped with ArcGIS 8 and appear not to have returned. Just for fun, here's a set of fast, simple calculations that will do the job no matter what.




      1. Convert the grid into a single zone by equating it with itself, as in



        "My grid" == "My grid"



      2. Create a column index grid by flow-accumulating a constant grid with value 1. (The indexes will start with 0.) If desired, multiply this by the cellsize and add the x-coordinate of the origin to obtain an x-coordinate grid "X" (shown below).


      3. Similarly, create a row index grid (and then a y-coordinate grid "Y") by flow-accumulating a constant grid with value 64.


      4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y": you now have the desired extent.


      Final image



      (The extent, as shown in the two tables of zonal statistics, is depicted by a rectangular outline in this figure. Grid "I" is the zone grid obtained in step (1).)



      To go further, you will need to extract these four numbers from their output tables and use them to limit the analysis extent. Copying the original grid, with the restricted analysis extent in place, completes the task.






      share|improve this answer



























        8














        The trick is to compute the limits of the data that have values. Perhaps the fastest, most natural, and most general way to obtain these is with zonal summaries: by using all non-NoData cells for the zone, the zonal min and max of grids containing the X and Y coordinates will provide the full extent.



        ESRI keeps changing the ways in which these calculations can be done; for instance, built-in expressions for the coordinate grids were dropped with ArcGIS 8 and appear not to have returned. Just for fun, here's a set of fast, simple calculations that will do the job no matter what.




        1. Convert the grid into a single zone by equating it with itself, as in



          "My grid" == "My grid"



        2. Create a column index grid by flow-accumulating a constant grid with value 1. (The indexes will start with 0.) If desired, multiply this by the cellsize and add the x-coordinate of the origin to obtain an x-coordinate grid "X" (shown below).


        3. Similarly, create a row index grid (and then a y-coordinate grid "Y") by flow-accumulating a constant grid with value 64.


        4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y": you now have the desired extent.


        Final image



        (The extent, as shown in the two tables of zonal statistics, is depicted by a rectangular outline in this figure. Grid "I" is the zone grid obtained in step (1).)



        To go further, you will need to extract these four numbers from their output tables and use them to limit the analysis extent. Copying the original grid, with the restricted analysis extent in place, completes the task.






        share|improve this answer

























          8












          8








          8







          The trick is to compute the limits of the data that have values. Perhaps the fastest, most natural, and most general way to obtain these is with zonal summaries: by using all non-NoData cells for the zone, the zonal min and max of grids containing the X and Y coordinates will provide the full extent.



          ESRI keeps changing the ways in which these calculations can be done; for instance, built-in expressions for the coordinate grids were dropped with ArcGIS 8 and appear not to have returned. Just for fun, here's a set of fast, simple calculations that will do the job no matter what.




          1. Convert the grid into a single zone by equating it with itself, as in



            "My grid" == "My grid"



          2. Create a column index grid by flow-accumulating a constant grid with value 1. (The indexes will start with 0.) If desired, multiply this by the cellsize and add the x-coordinate of the origin to obtain an x-coordinate grid "X" (shown below).


          3. Similarly, create a row index grid (and then a y-coordinate grid "Y") by flow-accumulating a constant grid with value 64.


          4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y": you now have the desired extent.


          Final image



          (The extent, as shown in the two tables of zonal statistics, is depicted by a rectangular outline in this figure. Grid "I" is the zone grid obtained in step (1).)



          To go further, you will need to extract these four numbers from their output tables and use them to limit the analysis extent. Copying the original grid, with the restricted analysis extent in place, completes the task.






          share|improve this answer













          The trick is to compute the limits of the data that have values. Perhaps the fastest, most natural, and most general way to obtain these is with zonal summaries: by using all non-NoData cells for the zone, the zonal min and max of grids containing the X and Y coordinates will provide the full extent.



          ESRI keeps changing the ways in which these calculations can be done; for instance, built-in expressions for the coordinate grids were dropped with ArcGIS 8 and appear not to have returned. Just for fun, here's a set of fast, simple calculations that will do the job no matter what.




          1. Convert the grid into a single zone by equating it with itself, as in



            "My grid" == "My grid"



          2. Create a column index grid by flow-accumulating a constant grid with value 1. (The indexes will start with 0.) If desired, multiply this by the cellsize and add the x-coordinate of the origin to obtain an x-coordinate grid "X" (shown below).


          3. Similarly, create a row index grid (and then a y-coordinate grid "Y") by flow-accumulating a constant grid with value 64.


          4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y": you now have the desired extent.


          Final image



          (The extent, as shown in the two tables of zonal statistics, is depicted by a rectangular outline in this figure. Grid "I" is the zone grid obtained in step (1).)



          To go further, you will need to extract these four numbers from their output tables and use them to limit the analysis extent. Copying the original grid, with the restricted analysis extent in place, completes the task.







          share|improve this answer












          share|improve this answer



          share|improve this answer










          answered Mar 21 '13 at 21:33









          whuberwhuber

          60.2k11156248




          60.2k11156248





















              8














              Here's a version of @whubers method for ArcGIS 10.1+ as a python toolbox (.pyt).



              import arcpy

              class Toolbox(object):
              def __init__(self):
              """Define the toolbox (the name of the toolbox is the name of the
              .pyt file)."""
              self.label = "Raster Toolbox"
              self.alias = ""

              # List of tool classes associated with this toolbox
              self.tools = [ClipNoData]


              class ClipNoData(object):
              def __init__(self):
              """Clip raster extent to the data that have values"""
              self.label = "Clip NoData"
              self.description = "Clip raster extent to the data that have values. "
              self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

              self.canRunInBackground = False

              def getParameterInfo(self):
              """Define parameter definitions"""
              params = []

              # First parameter
              params+=[arcpy.Parameter(
              displayName="Input Raster",
              name="in_raster",
              datatype='GPRasterLayer',
              parameterType="Required",
              direction="Input")
              ]

              # Second parameter
              params+=[arcpy.Parameter(
              displayName="Output Raster",
              name="out_raster",
              datatype="DERasterDataset",
              parameterType="Required",
              direction="Output")
              ]

              return params

              def isLicensed(self):
              """Set whether tool is licensed to execute."""
              return arcpy.CheckExtension('spatial')==u'Available'

              def execute(self, parameters, messages):
              """See https://gis.stackexchange.com/a/55150/2856
              ##Code comments paraphrased from @whubers GIS StackExchange answer
              """
              try:
              #Setup
              arcpy.CheckOutExtension('spatial')
              from arcpy.sa import *
              in_raster = parameters[0].valueAsText
              out_raster = parameters[1].valueAsText

              dsc=arcpy.Describe(in_raster)
              xmin=dsc.extent.XMin
              ymin=dsc.extent.YMin
              mx=dsc.meanCellWidth
              my=dsc.meanCellHeight
              arcpy.env.extent=in_raster
              arcpy.env.cellSize=in_raster
              arcpy.AddMessage(out_raster)

              ## 1. Convert the grid into a single zone by equating it with itself
              arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
              zones = Raster(in_raster) == Raster(in_raster)

              ## 2. Create a column index grid by flow-accumulating a constant grid
              ## with value 1. (The indexes will start with 0.) Multiply this by
              ## the cellsize and add the x-coordinate of the origin to obtain
              ## an x-coordinate grid.
              arcpy.AddMessage(r'Create a constant grid...')
              const = CreateConstantRaster(1)

              arcpy.AddMessage(r'2. Create an x-coordinate grid...')
              xmap = (FlowAccumulation(const)) * mx + xmin

              ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
              ## constant grid with value 64.
              arcpy.AddMessage(r'3. Create a y-coordinate grid...')
              ymap = (FlowAccumulation(const * 64)) * my + ymin

              ## 4. Use the zone grid from step (1) to compute the zonal min and
              ## max of "X" and "Y"
              arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

              xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORYxrange", "NODATA", "MIN_MAX")
              xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORYxrange", ["MIN","MAX"]).next()

              yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORYyrange", "NODATA", "MIN_MAX")
              ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORYyrange", ["MIN","MAX"]).next()

              arcpy.Delete_management(r"IN_MEMORYxrange")
              arcpy.Delete_management(r"IN_MEMORYyrange")

              # Write out the reduced raster
              arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
              output = Raster(in_raster) * 1
              output.save(out_raster)

              except:raise
              finally:arcpy.CheckInExtension('spatial')





              share|improve this answer

























              • Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

                – matt wilkie
                Mar 22 '13 at 15:52






              • 1





                I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

                – user2856
                Mar 22 '13 at 20:31











              • good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

                – matt wilkie
                Mar 22 '13 at 20:42











              • +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

                – whuber
                Mar 23 '13 at 3:34















              8














              Here's a version of @whubers method for ArcGIS 10.1+ as a python toolbox (.pyt).



              import arcpy

              class Toolbox(object):
              def __init__(self):
              """Define the toolbox (the name of the toolbox is the name of the
              .pyt file)."""
              self.label = "Raster Toolbox"
              self.alias = ""

              # List of tool classes associated with this toolbox
              self.tools = [ClipNoData]


              class ClipNoData(object):
              def __init__(self):
              """Clip raster extent to the data that have values"""
              self.label = "Clip NoData"
              self.description = "Clip raster extent to the data that have values. "
              self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

              self.canRunInBackground = False

              def getParameterInfo(self):
              """Define parameter definitions"""
              params = []

              # First parameter
              params+=[arcpy.Parameter(
              displayName="Input Raster",
              name="in_raster",
              datatype='GPRasterLayer',
              parameterType="Required",
              direction="Input")
              ]

              # Second parameter
              params+=[arcpy.Parameter(
              displayName="Output Raster",
              name="out_raster",
              datatype="DERasterDataset",
              parameterType="Required",
              direction="Output")
              ]

              return params

              def isLicensed(self):
              """Set whether tool is licensed to execute."""
              return arcpy.CheckExtension('spatial')==u'Available'

              def execute(self, parameters, messages):
              """See https://gis.stackexchange.com/a/55150/2856
              ##Code comments paraphrased from @whubers GIS StackExchange answer
              """
              try:
              #Setup
              arcpy.CheckOutExtension('spatial')
              from arcpy.sa import *
              in_raster = parameters[0].valueAsText
              out_raster = parameters[1].valueAsText

              dsc=arcpy.Describe(in_raster)
              xmin=dsc.extent.XMin
              ymin=dsc.extent.YMin
              mx=dsc.meanCellWidth
              my=dsc.meanCellHeight
              arcpy.env.extent=in_raster
              arcpy.env.cellSize=in_raster
              arcpy.AddMessage(out_raster)

              ## 1. Convert the grid into a single zone by equating it with itself
              arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
              zones = Raster(in_raster) == Raster(in_raster)

              ## 2. Create a column index grid by flow-accumulating a constant grid
              ## with value 1. (The indexes will start with 0.) Multiply this by
              ## the cellsize and add the x-coordinate of the origin to obtain
              ## an x-coordinate grid.
              arcpy.AddMessage(r'Create a constant grid...')
              const = CreateConstantRaster(1)

              arcpy.AddMessage(r'2. Create an x-coordinate grid...')
              xmap = (FlowAccumulation(const)) * mx + xmin

              ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
              ## constant grid with value 64.
              arcpy.AddMessage(r'3. Create a y-coordinate grid...')
              ymap = (FlowAccumulation(const * 64)) * my + ymin

              ## 4. Use the zone grid from step (1) to compute the zonal min and
              ## max of "X" and "Y"
              arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

              xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORYxrange", "NODATA", "MIN_MAX")
              xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORYxrange", ["MIN","MAX"]).next()

              yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORYyrange", "NODATA", "MIN_MAX")
              ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORYyrange", ["MIN","MAX"]).next()

              arcpy.Delete_management(r"IN_MEMORYxrange")
              arcpy.Delete_management(r"IN_MEMORYyrange")

              # Write out the reduced raster
              arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
              output = Raster(in_raster) * 1
              output.save(out_raster)

              except:raise
              finally:arcpy.CheckInExtension('spatial')





              share|improve this answer

























              • Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

                – matt wilkie
                Mar 22 '13 at 15:52






              • 1





                I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

                – user2856
                Mar 22 '13 at 20:31











              • good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

                – matt wilkie
                Mar 22 '13 at 20:42











              • +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

                – whuber
                Mar 23 '13 at 3:34













              8












              8








              8







              Here's a version of @whubers method for ArcGIS 10.1+ as a python toolbox (.pyt).



              import arcpy

              class Toolbox(object):
              def __init__(self):
              """Define the toolbox (the name of the toolbox is the name of the
              .pyt file)."""
              self.label = "Raster Toolbox"
              self.alias = ""

              # List of tool classes associated with this toolbox
              self.tools = [ClipNoData]


              class ClipNoData(object):
              def __init__(self):
              """Clip raster extent to the data that have values"""
              self.label = "Clip NoData"
              self.description = "Clip raster extent to the data that have values. "
              self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

              self.canRunInBackground = False

              def getParameterInfo(self):
              """Define parameter definitions"""
              params = []

              # First parameter
              params+=[arcpy.Parameter(
              displayName="Input Raster",
              name="in_raster",
              datatype='GPRasterLayer',
              parameterType="Required",
              direction="Input")
              ]

              # Second parameter
              params+=[arcpy.Parameter(
              displayName="Output Raster",
              name="out_raster",
              datatype="DERasterDataset",
              parameterType="Required",
              direction="Output")
              ]

              return params

              def isLicensed(self):
              """Set whether tool is licensed to execute."""
              return arcpy.CheckExtension('spatial')==u'Available'

              def execute(self, parameters, messages):
              """See https://gis.stackexchange.com/a/55150/2856
              ##Code comments paraphrased from @whubers GIS StackExchange answer
              """
              try:
              #Setup
              arcpy.CheckOutExtension('spatial')
              from arcpy.sa import *
              in_raster = parameters[0].valueAsText
              out_raster = parameters[1].valueAsText

              dsc=arcpy.Describe(in_raster)
              xmin=dsc.extent.XMin
              ymin=dsc.extent.YMin
              mx=dsc.meanCellWidth
              my=dsc.meanCellHeight
              arcpy.env.extent=in_raster
              arcpy.env.cellSize=in_raster
              arcpy.AddMessage(out_raster)

              ## 1. Convert the grid into a single zone by equating it with itself
              arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
              zones = Raster(in_raster) == Raster(in_raster)

              ## 2. Create a column index grid by flow-accumulating a constant grid
              ## with value 1. (The indexes will start with 0.) Multiply this by
              ## the cellsize and add the x-coordinate of the origin to obtain
              ## an x-coordinate grid.
              arcpy.AddMessage(r'Create a constant grid...')
              const = CreateConstantRaster(1)

              arcpy.AddMessage(r'2. Create an x-coordinate grid...')
              xmap = (FlowAccumulation(const)) * mx + xmin

              ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
              ## constant grid with value 64.
              arcpy.AddMessage(r'3. Create a y-coordinate grid...')
              ymap = (FlowAccumulation(const * 64)) * my + ymin

              ## 4. Use the zone grid from step (1) to compute the zonal min and
              ## max of "X" and "Y"
              arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

              xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORYxrange", "NODATA", "MIN_MAX")
              xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORYxrange", ["MIN","MAX"]).next()

              yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORYyrange", "NODATA", "MIN_MAX")
              ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORYyrange", ["MIN","MAX"]).next()

              arcpy.Delete_management(r"IN_MEMORYxrange")
              arcpy.Delete_management(r"IN_MEMORYyrange")

              # Write out the reduced raster
              arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
              output = Raster(in_raster) * 1
              output.save(out_raster)

              except:raise
              finally:arcpy.CheckInExtension('spatial')





              share|improve this answer















              Here's a version of @whubers method for ArcGIS 10.1+ as a python toolbox (.pyt).



              import arcpy

              class Toolbox(object):
              def __init__(self):
              """Define the toolbox (the name of the toolbox is the name of the
              .pyt file)."""
              self.label = "Raster Toolbox"
              self.alias = ""

              # List of tool classes associated with this toolbox
              self.tools = [ClipNoData]


              class ClipNoData(object):
              def __init__(self):
              """Clip raster extent to the data that have values"""
              self.label = "Clip NoData"
              self.description = "Clip raster extent to the data that have values. "
              self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

              self.canRunInBackground = False

              def getParameterInfo(self):
              """Define parameter definitions"""
              params = []

              # First parameter
              params+=[arcpy.Parameter(
              displayName="Input Raster",
              name="in_raster",
              datatype='GPRasterLayer',
              parameterType="Required",
              direction="Input")
              ]

              # Second parameter
              params+=[arcpy.Parameter(
              displayName="Output Raster",
              name="out_raster",
              datatype="DERasterDataset",
              parameterType="Required",
              direction="Output")
              ]

              return params

              def isLicensed(self):
              """Set whether tool is licensed to execute."""
              return arcpy.CheckExtension('spatial')==u'Available'

              def execute(self, parameters, messages):
              """See https://gis.stackexchange.com/a/55150/2856
              ##Code comments paraphrased from @whubers GIS StackExchange answer
              """
              try:
              #Setup
              arcpy.CheckOutExtension('spatial')
              from arcpy.sa import *
              in_raster = parameters[0].valueAsText
              out_raster = parameters[1].valueAsText

              dsc=arcpy.Describe(in_raster)
              xmin=dsc.extent.XMin
              ymin=dsc.extent.YMin
              mx=dsc.meanCellWidth
              my=dsc.meanCellHeight
              arcpy.env.extent=in_raster
              arcpy.env.cellSize=in_raster
              arcpy.AddMessage(out_raster)

              ## 1. Convert the grid into a single zone by equating it with itself
              arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
              zones = Raster(in_raster) == Raster(in_raster)

              ## 2. Create a column index grid by flow-accumulating a constant grid
              ## with value 1. (The indexes will start with 0.) Multiply this by
              ## the cellsize and add the x-coordinate of the origin to obtain
              ## an x-coordinate grid.
              arcpy.AddMessage(r'Create a constant grid...')
              const = CreateConstantRaster(1)

              arcpy.AddMessage(r'2. Create an x-coordinate grid...')
              xmap = (FlowAccumulation(const)) * mx + xmin

              ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
              ## constant grid with value 64.
              arcpy.AddMessage(r'3. Create a y-coordinate grid...')
              ymap = (FlowAccumulation(const * 64)) * my + ymin

              ## 4. Use the zone grid from step (1) to compute the zonal min and
              ## max of "X" and "Y"
              arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

              xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORYxrange", "NODATA", "MIN_MAX")
              xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORYxrange", ["MIN","MAX"]).next()

              yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORYyrange", "NODATA", "MIN_MAX")
              ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORYyrange", ["MIN","MAX"]).next()

              arcpy.Delete_management(r"IN_MEMORYxrange")
              arcpy.Delete_management(r"IN_MEMORYyrange")

              # Write out the reduced raster
              arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
              output = Raster(in_raster) * 1
              output.save(out_raster)

              except:raise
              finally:arcpy.CheckInExtension('spatial')






              share|improve this answer














              share|improve this answer



              share|improve this answer








              edited Apr 13 '17 at 12:33









              Community

              1




              1










              answered Mar 22 '13 at 4:44









              user2856user2856

              30.6k258106




              30.6k258106












              • Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

                – matt wilkie
                Mar 22 '13 at 15:52






              • 1





                I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

                – user2856
                Mar 22 '13 at 20:31











              • good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

                – matt wilkie
                Mar 22 '13 at 20:42











              • +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

                – whuber
                Mar 23 '13 at 3:34

















              • Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

                – matt wilkie
                Mar 22 '13 at 15:52






              • 1





                I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

                – user2856
                Mar 22 '13 at 20:31











              • good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

                – matt wilkie
                Mar 22 '13 at 20:42











              • +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

                – whuber
                Mar 23 '13 at 3:34
















              Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

              – matt wilkie
              Mar 22 '13 at 15:52





              Very nice Luke. Self contained, runnable, uses in_memory, and well commented to boot. I had to turn off background processing (Geoprocessing > options >...) to get it working.

              – matt wilkie
              Mar 22 '13 at 15:52




              1




              1





              I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

              – user2856
              Mar 22 '13 at 20:31





              I've updated the script and set canRunInBackground=False. I will note that it's worth changing the workspace/scratchworkspace environments to a local folder (not FGDB) as I found when I left them as default (i.e <network user profile>Default.gdb) the script took 9min!!! to run on a 250x250 cell grid. Changing to a local FGDB it took 9sec and to a local folder 1-2sec...

              – user2856
              Mar 22 '13 at 20:31













              good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

              – matt wilkie
              Mar 22 '13 at 20:42





              good point about local folders, and thanks for the quick background fix (much better than writing instructions for everyone I pass it to). Might be worth throwing this up on bitbucket/github/gcode/etc.

              – matt wilkie
              Mar 22 '13 at 20:42













              +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

              – whuber
              Mar 23 '13 at 3:34





              +1 Thanks for this contribution, Luke! I appreciate you filling in the (rather large) gap left in my answer.

              – whuber
              Mar 23 '13 at 3:34











              4














              I've devised a gdal and numpy based solution. It breaks the raster matrix into rows and columns and drop any empty row/column. In this implementation "empty" is anything less than 1, and only single band rasters accounted for.



              (I realise as I write that this scanline approach is only suitable for images with nodata "collars". If your data are islands in seas of nulls, the space between islands will get dropped as well, squishing everything together and totally messing up the georeferencing.)



              The business parts (needs fleshing out, won't work as is):





               #read raster into a numpy array
              data = np.array(gdal.Open(src_raster).ReadAsArray())
              #scan for data
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              # assumes data is any value greater than zero
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # retrieve source geo reference info
              georef = raster.GetGeoTransform()
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              # write to disk
              band = out_raster.GetRasterBand(1)
              band.WriteArray(new_data)
              band.FlushCache()
              out_raster = None


              In a full script:



              import os
              import sys
              import numpy as np
              from osgeo import gdal

              if len(sys.argv) < 2:
              print 'n [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
              sys.exit(1)

              src_raster = sys.argv[1]
              out_raster = sys.argv[2]

              def main(src_raster):
              raster = gdal.Open(src_raster)

              # Read georeferencing, oriented from top-left
              # ref:GDAL Tutorial, Getting Dataset Information
              georef = raster.GetGeoTransform()
              print 'nSource raster (geo units):'
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]
              cols, rows = raster.RasterYSize, raster.RasterXSize
              print ' Origin (top left): :10, :10'.format(xmin, ymax)
              print ' Pixel size (x,-y): :10, :10'.format(xcell, ycell)
              print ' Columns, rows : :10, :10'.format(cols, rows)

              # Transfer to numpy and scan for data
              # oriented from bottom-left
              data = np.array(raster.ReadAsArray())
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              new_rows, new_cols = new_data.shape # note: inverted relative to geo units
              #print cropped_transform

              print 'nCrop box (pixel units):', crop_box
              print ' Stripped columns : :10'.format(cols - new_cols)
              print ' Stripped rows : :10'.format(rows - new_rows)

              print 'nCropped raster (geo units):'
              print ' Origin (top left): :10, :10'.format(new_xmin, new_ymax)
              print ' Columns, rows : :10, :10'.format(new_cols, new_rows)

              raster = None
              return new_data, cropped_transform


              def write_raster(template, array, transform, filename):
              '''Create a new raster from an array.

              template = raster dataset to copy projection info from
              array = numpy array of a raster
              transform = geo referencing (x,y origin and pixel dimensions)
              filename = path to output image (will be overwritten)
              '''
              template = gdal.Open(template)
              driver = template.GetDriver()
              rows,cols = array.shape
              out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
              out_raster.SetGeoTransform(transform)
              out_raster.SetProjection(template.GetProjection())
              band = out_raster.GetRasterBand(1)
              band.WriteArray(array)
              band.FlushCache()
              out_raster = None
              template = None

              if __name__ == '__main__':
              cropped_raster, cropped_transform = main(src_raster)
              write_raster(src_raster, cropped_raster, cropped_transform, out_raster)


              The script is in my code stash on Github, if the link goes 404 hunt around a bit; these folders are ripe for some reorganization.






              share|improve this answer




















              • 1





                This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

                – user32882
                Oct 24 '18 at 13:05
















              4














              I've devised a gdal and numpy based solution. It breaks the raster matrix into rows and columns and drop any empty row/column. In this implementation "empty" is anything less than 1, and only single band rasters accounted for.



              (I realise as I write that this scanline approach is only suitable for images with nodata "collars". If your data are islands in seas of nulls, the space between islands will get dropped as well, squishing everything together and totally messing up the georeferencing.)



              The business parts (needs fleshing out, won't work as is):





               #read raster into a numpy array
              data = np.array(gdal.Open(src_raster).ReadAsArray())
              #scan for data
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              # assumes data is any value greater than zero
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # retrieve source geo reference info
              georef = raster.GetGeoTransform()
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              # write to disk
              band = out_raster.GetRasterBand(1)
              band.WriteArray(new_data)
              band.FlushCache()
              out_raster = None


              In a full script:



              import os
              import sys
              import numpy as np
              from osgeo import gdal

              if len(sys.argv) < 2:
              print 'n [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
              sys.exit(1)

              src_raster = sys.argv[1]
              out_raster = sys.argv[2]

              def main(src_raster):
              raster = gdal.Open(src_raster)

              # Read georeferencing, oriented from top-left
              # ref:GDAL Tutorial, Getting Dataset Information
              georef = raster.GetGeoTransform()
              print 'nSource raster (geo units):'
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]
              cols, rows = raster.RasterYSize, raster.RasterXSize
              print ' Origin (top left): :10, :10'.format(xmin, ymax)
              print ' Pixel size (x,-y): :10, :10'.format(xcell, ycell)
              print ' Columns, rows : :10, :10'.format(cols, rows)

              # Transfer to numpy and scan for data
              # oriented from bottom-left
              data = np.array(raster.ReadAsArray())
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              new_rows, new_cols = new_data.shape # note: inverted relative to geo units
              #print cropped_transform

              print 'nCrop box (pixel units):', crop_box
              print ' Stripped columns : :10'.format(cols - new_cols)
              print ' Stripped rows : :10'.format(rows - new_rows)

              print 'nCropped raster (geo units):'
              print ' Origin (top left): :10, :10'.format(new_xmin, new_ymax)
              print ' Columns, rows : :10, :10'.format(new_cols, new_rows)

              raster = None
              return new_data, cropped_transform


              def write_raster(template, array, transform, filename):
              '''Create a new raster from an array.

              template = raster dataset to copy projection info from
              array = numpy array of a raster
              transform = geo referencing (x,y origin and pixel dimensions)
              filename = path to output image (will be overwritten)
              '''
              template = gdal.Open(template)
              driver = template.GetDriver()
              rows,cols = array.shape
              out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
              out_raster.SetGeoTransform(transform)
              out_raster.SetProjection(template.GetProjection())
              band = out_raster.GetRasterBand(1)
              band.WriteArray(array)
              band.FlushCache()
              out_raster = None
              template = None

              if __name__ == '__main__':
              cropped_raster, cropped_transform = main(src_raster)
              write_raster(src_raster, cropped_raster, cropped_transform, out_raster)


              The script is in my code stash on Github, if the link goes 404 hunt around a bit; these folders are ripe for some reorganization.






              share|improve this answer




















              • 1





                This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

                – user32882
                Oct 24 '18 at 13:05














              4












              4








              4







              I've devised a gdal and numpy based solution. It breaks the raster matrix into rows and columns and drop any empty row/column. In this implementation "empty" is anything less than 1, and only single band rasters accounted for.



              (I realise as I write that this scanline approach is only suitable for images with nodata "collars". If your data are islands in seas of nulls, the space between islands will get dropped as well, squishing everything together and totally messing up the georeferencing.)



              The business parts (needs fleshing out, won't work as is):





               #read raster into a numpy array
              data = np.array(gdal.Open(src_raster).ReadAsArray())
              #scan for data
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              # assumes data is any value greater than zero
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # retrieve source geo reference info
              georef = raster.GetGeoTransform()
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              # write to disk
              band = out_raster.GetRasterBand(1)
              band.WriteArray(new_data)
              band.FlushCache()
              out_raster = None


              In a full script:



              import os
              import sys
              import numpy as np
              from osgeo import gdal

              if len(sys.argv) < 2:
              print 'n [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
              sys.exit(1)

              src_raster = sys.argv[1]
              out_raster = sys.argv[2]

              def main(src_raster):
              raster = gdal.Open(src_raster)

              # Read georeferencing, oriented from top-left
              # ref:GDAL Tutorial, Getting Dataset Information
              georef = raster.GetGeoTransform()
              print 'nSource raster (geo units):'
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]
              cols, rows = raster.RasterYSize, raster.RasterXSize
              print ' Origin (top left): :10, :10'.format(xmin, ymax)
              print ' Pixel size (x,-y): :10, :10'.format(xcell, ycell)
              print ' Columns, rows : :10, :10'.format(cols, rows)

              # Transfer to numpy and scan for data
              # oriented from bottom-left
              data = np.array(raster.ReadAsArray())
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              new_rows, new_cols = new_data.shape # note: inverted relative to geo units
              #print cropped_transform

              print 'nCrop box (pixel units):', crop_box
              print ' Stripped columns : :10'.format(cols - new_cols)
              print ' Stripped rows : :10'.format(rows - new_rows)

              print 'nCropped raster (geo units):'
              print ' Origin (top left): :10, :10'.format(new_xmin, new_ymax)
              print ' Columns, rows : :10, :10'.format(new_cols, new_rows)

              raster = None
              return new_data, cropped_transform


              def write_raster(template, array, transform, filename):
              '''Create a new raster from an array.

              template = raster dataset to copy projection info from
              array = numpy array of a raster
              transform = geo referencing (x,y origin and pixel dimensions)
              filename = path to output image (will be overwritten)
              '''
              template = gdal.Open(template)
              driver = template.GetDriver()
              rows,cols = array.shape
              out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
              out_raster.SetGeoTransform(transform)
              out_raster.SetProjection(template.GetProjection())
              band = out_raster.GetRasterBand(1)
              band.WriteArray(array)
              band.FlushCache()
              out_raster = None
              template = None

              if __name__ == '__main__':
              cropped_raster, cropped_transform = main(src_raster)
              write_raster(src_raster, cropped_raster, cropped_transform, out_raster)


              The script is in my code stash on Github, if the link goes 404 hunt around a bit; these folders are ripe for some reorganization.






              share|improve this answer















              I've devised a gdal and numpy based solution. It breaks the raster matrix into rows and columns and drop any empty row/column. In this implementation "empty" is anything less than 1, and only single band rasters accounted for.



              (I realise as I write that this scanline approach is only suitable for images with nodata "collars". If your data are islands in seas of nulls, the space between islands will get dropped as well, squishing everything together and totally messing up the georeferencing.)



              The business parts (needs fleshing out, won't work as is):





               #read raster into a numpy array
              data = np.array(gdal.Open(src_raster).ReadAsArray())
              #scan for data
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              # assumes data is any value greater than zero
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # retrieve source geo reference info
              georef = raster.GetGeoTransform()
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              # write to disk
              band = out_raster.GetRasterBand(1)
              band.WriteArray(new_data)
              band.FlushCache()
              out_raster = None


              In a full script:



              import os
              import sys
              import numpy as np
              from osgeo import gdal

              if len(sys.argv) < 2:
              print 'n [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
              sys.exit(1)

              src_raster = sys.argv[1]
              out_raster = sys.argv[2]

              def main(src_raster):
              raster = gdal.Open(src_raster)

              # Read georeferencing, oriented from top-left
              # ref:GDAL Tutorial, Getting Dataset Information
              georef = raster.GetGeoTransform()
              print 'nSource raster (geo units):'
              xmin, ymax = georef[0], georef[3]
              xcell, ycell = georef[1], georef[5]
              cols, rows = raster.RasterYSize, raster.RasterXSize
              print ' Origin (top left): :10, :10'.format(xmin, ymax)
              print ' Pixel size (x,-y): :10, :10'.format(xcell, ycell)
              print ' Columns, rows : :10, :10'.format(cols, rows)

              # Transfer to numpy and scan for data
              # oriented from bottom-left
              data = np.array(raster.ReadAsArray())
              non_empty_columns = np.where(data.max(axis=0)>0)[0]
              non_empty_rows = np.where(data.max(axis=1)>0)[0]
              crop_box = (min(non_empty_rows), max(non_empty_rows),
              min(non_empty_columns), max(non_empty_columns))

              # Calculate cropped geo referencing
              new_xmin = xmin + (xcell * crop_box[0]) + xcell
              new_ymax = ymax + (ycell * crop_box[2]) - ycell
              cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

              # crop
              new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

              new_rows, new_cols = new_data.shape # note: inverted relative to geo units
              #print cropped_transform

              print 'nCrop box (pixel units):', crop_box
              print ' Stripped columns : :10'.format(cols - new_cols)
              print ' Stripped rows : :10'.format(rows - new_rows)

              print 'nCropped raster (geo units):'
              print ' Origin (top left): :10, :10'.format(new_xmin, new_ymax)
              print ' Columns, rows : :10, :10'.format(new_cols, new_rows)

              raster = None
              return new_data, cropped_transform


              def write_raster(template, array, transform, filename):
              '''Create a new raster from an array.

              template = raster dataset to copy projection info from
              array = numpy array of a raster
              transform = geo referencing (x,y origin and pixel dimensions)
              filename = path to output image (will be overwritten)
              '''
              template = gdal.Open(template)
              driver = template.GetDriver()
              rows,cols = array.shape
              out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
              out_raster.SetGeoTransform(transform)
              out_raster.SetProjection(template.GetProjection())
              band = out_raster.GetRasterBand(1)
              band.WriteArray(array)
              band.FlushCache()
              out_raster = None
              template = None

              if __name__ == '__main__':
              cropped_raster, cropped_transform = main(src_raster)
              write_raster(src_raster, cropped_raster, cropped_transform, out_raster)


              The script is in my code stash on Github, if the link goes 404 hunt around a bit; these folders are ripe for some reorganization.







              share|improve this answer














              share|improve this answer



              share|improve this answer








              edited yesterday









              Martin Valgur

              1,0991015




              1,0991015










              answered Jun 2 '15 at 18:51









              matt wilkiematt wilkie

              18.2k30116241




              18.2k30116241







              • 1





                This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

                – user32882
                Oct 24 '18 at 13:05













              • 1





                This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

                – user32882
                Oct 24 '18 at 13:05








              1




              1





              This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

              – user32882
              Oct 24 '18 at 13:05






              This won'y work for really large datasets. I get a MemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError

              – user32882
              Oct 24 '18 at 13:05












              2














              For all its analytical power, ArcGIS lacks basic raster manipulations that you can find with traditional desktop image editors like GIMP. It expects that you want to use the same analysis extent for your output raster as your input raster, unless you manually override the Output Extent environment setting. Since this is exactly what you are looking to find, not set, the ArcGIS way of doing things is getting in the way.



              Despite my best efforts, and without resorting to programming, I could find no way to get the extent of your desired subset of the image (without raster-to-vector conversion which is computationally wasteful).



              Instead, I used GIMP to select the blue area using the select by color tool and then inverted the selection, hit Delete to remove the rest of the pixels, inverted the selection again, cropped the image to selection, and finally exported it back out to PNG. GIMP saved it as a 1-bit depth image. The result is below:



              Output



              Of course, because your sample image lacked a spatial reference component, and GIMP is not spatially aware, the output image is about as useful as your sample input. You will need to georeference it for it to be of use in a spatial analysis.






              share|improve this answer























              • Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

                – whuber
                Jan 4 '13 at 14:45











              • @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

                – Seen
                Jan 4 '13 at 15:24











              • @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

                – whuber
                Jan 4 '13 at 16:00















              2














              For all its analytical power, ArcGIS lacks basic raster manipulations that you can find with traditional desktop image editors like GIMP. It expects that you want to use the same analysis extent for your output raster as your input raster, unless you manually override the Output Extent environment setting. Since this is exactly what you are looking to find, not set, the ArcGIS way of doing things is getting in the way.



              Despite my best efforts, and without resorting to programming, I could find no way to get the extent of your desired subset of the image (without raster-to-vector conversion which is computationally wasteful).



              Instead, I used GIMP to select the blue area using the select by color tool and then inverted the selection, hit Delete to remove the rest of the pixels, inverted the selection again, cropped the image to selection, and finally exported it back out to PNG. GIMP saved it as a 1-bit depth image. The result is below:



              Output



              Of course, because your sample image lacked a spatial reference component, and GIMP is not spatially aware, the output image is about as useful as your sample input. You will need to georeference it for it to be of use in a spatial analysis.






              share|improve this answer























              • Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

                – whuber
                Jan 4 '13 at 14:45











              • @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

                – Seen
                Jan 4 '13 at 15:24











              • @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

                – whuber
                Jan 4 '13 at 16:00













              2












              2








              2







              For all its analytical power, ArcGIS lacks basic raster manipulations that you can find with traditional desktop image editors like GIMP. It expects that you want to use the same analysis extent for your output raster as your input raster, unless you manually override the Output Extent environment setting. Since this is exactly what you are looking to find, not set, the ArcGIS way of doing things is getting in the way.



              Despite my best efforts, and without resorting to programming, I could find no way to get the extent of your desired subset of the image (without raster-to-vector conversion which is computationally wasteful).



              Instead, I used GIMP to select the blue area using the select by color tool and then inverted the selection, hit Delete to remove the rest of the pixels, inverted the selection again, cropped the image to selection, and finally exported it back out to PNG. GIMP saved it as a 1-bit depth image. The result is below:



              Output



              Of course, because your sample image lacked a spatial reference component, and GIMP is not spatially aware, the output image is about as useful as your sample input. You will need to georeference it for it to be of use in a spatial analysis.






              share|improve this answer













              For all its analytical power, ArcGIS lacks basic raster manipulations that you can find with traditional desktop image editors like GIMP. It expects that you want to use the same analysis extent for your output raster as your input raster, unless you manually override the Output Extent environment setting. Since this is exactly what you are looking to find, not set, the ArcGIS way of doing things is getting in the way.



              Despite my best efforts, and without resorting to programming, I could find no way to get the extent of your desired subset of the image (without raster-to-vector conversion which is computationally wasteful).



              Instead, I used GIMP to select the blue area using the select by color tool and then inverted the selection, hit Delete to remove the rest of the pixels, inverted the selection again, cropped the image to selection, and finally exported it back out to PNG. GIMP saved it as a 1-bit depth image. The result is below:



              Output



              Of course, because your sample image lacked a spatial reference component, and GIMP is not spatially aware, the output image is about as useful as your sample input. You will need to georeference it for it to be of use in a spatial analysis.







              share|improve this answer












              share|improve this answer



              share|improve this answer










              answered Jan 4 '13 at 9:15









              blah238blah238

              32.3k677168




              32.3k677168












              • Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

                – whuber
                Jan 4 '13 at 14:45











              • @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

                – Seen
                Jan 4 '13 at 15:24











              • @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

                – whuber
                Jan 4 '13 at 16:00

















              • Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

                – whuber
                Jan 4 '13 at 14:45











              • @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

                – Seen
                Jan 4 '13 at 15:24











              • @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

                – whuber
                Jan 4 '13 at 16:00
















              Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

              – whuber
              Jan 4 '13 at 14:45





              Actually, this operation used to be easy in previous versions of Spatial Analyst: the zonal max and min of the two coordinate grids (X and Y), using the feature as the zone, give the extent exactly. (Well, you might want to expand it by half a cellsize in all four directions.) In ArcGIS 10, you need to be creative (or to use Python) to make a coordinate grid. Regardless, the whole thing can be done entirely within SA using only grid operations and no manual intervention.

              – whuber
              Jan 4 '13 at 14:45













              @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

              – Seen
              Jan 4 '13 at 15:24





              @whuber I saw your solution somewhere else, but still not sure how I can implement your method. Could you show me some more detail of this?

              – Seen
              Jan 4 '13 at 15:24













              @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

              – whuber
              Jan 4 '13 at 16:00





              @Seen A quick search of this site finds an account of the method at gis.stackexchange.com/a/13467.

              – whuber
              Jan 4 '13 at 16:00











              1














              Here is one possibility using SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021




              In SAGA GIS there is a "Crop to Data" module (in Grid Tools module
              library), which performs the task.



              But this would require you to import your Geotif with the GDAL import
              module, process it in SAGA, and finally export it as Geotif again with
              the GDAL export module.




              Another possibility using only ArcGIS GP tools would be to build a TIN from your raster using Raster to TIN, compute its boundary using TIN Domain, and Clip your raster by the boundary (or its envelope using Feature Envelope to Polygon).






              share|improve this answer





























                1














                Here is one possibility using SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021




                In SAGA GIS there is a "Crop to Data" module (in Grid Tools module
                library), which performs the task.



                But this would require you to import your Geotif with the GDAL import
                module, process it in SAGA, and finally export it as Geotif again with
                the GDAL export module.




                Another possibility using only ArcGIS GP tools would be to build a TIN from your raster using Raster to TIN, compute its boundary using TIN Domain, and Clip your raster by the boundary (or its envelope using Feature Envelope to Polygon).






                share|improve this answer



























                  1












                  1








                  1







                  Here is one possibility using SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021




                  In SAGA GIS there is a "Crop to Data" module (in Grid Tools module
                  library), which performs the task.



                  But this would require you to import your Geotif with the GDAL import
                  module, process it in SAGA, and finally export it as Geotif again with
                  the GDAL export module.




                  Another possibility using only ArcGIS GP tools would be to build a TIN from your raster using Raster to TIN, compute its boundary using TIN Domain, and Clip your raster by the boundary (or its envelope using Feature Envelope to Polygon).






                  share|improve this answer















                  Here is one possibility using SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021




                  In SAGA GIS there is a "Crop to Data" module (in Grid Tools module
                  library), which performs the task.



                  But this would require you to import your Geotif with the GDAL import
                  module, process it in SAGA, and finally export it as Geotif again with
                  the GDAL export module.




                  Another possibility using only ArcGIS GP tools would be to build a TIN from your raster using Raster to TIN, compute its boundary using TIN Domain, and Clip your raster by the boundary (or its envelope using Feature Envelope to Polygon).







                  share|improve this answer














                  share|improve this answer



                  share|improve this answer








                  edited Mar 21 '13 at 21:41

























                  answered Mar 21 '13 at 21:28









                  blah238blah238

                  32.3k677168




                  32.3k677168



























                      draft saved

                      draft discarded
















































                      Thanks for contributing an answer to Geographic Information Systems Stack Exchange!


                      • Please be sure to answer the question. Provide details and share your research!

                      But avoid


                      • Asking for help, clarification, or responding to other answers.

                      • Making statements based on opinion; back them up with references or personal experience.

                      To learn more, see our tips on writing great answers.




                      draft saved


                      draft discarded














                      StackExchange.ready(
                      function ()
                      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fgis.stackexchange.com%2fquestions%2f45159%2ffinding-minimum-bounding-extent-of-given-pixel-value-within-raster%23new-answer', 'question_page');

                      );

                      Post as a guest















                      Required, but never shown





















































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown

































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown







                      Popular posts from this blog

                      រឿង រ៉ូមេអូ និង ហ្ស៊ុយលីយេ សង្ខេបរឿង តួអង្គ បញ្ជីណែនាំ

                      Crop image to path created in TikZ? Announcing the arrival of Valued Associate #679: Cesar Manara Planned maintenance scheduled April 17/18, 2019 at 00:00UTC (8:00pm US/Eastern)Crop an inserted image?TikZ pictures does not appear in posterImage behind and beyond crop marks?Tikz picture as large as possible on A4 PageTransparency vs image compression dilemmaHow to crop background from image automatically?Image does not cropTikzexternal capturing crop marks when externalizing pgfplots?How to include image path that contains a dollar signCrop image with left size given

                      Romeo and Juliet ContentsCharactersSynopsisSourcesDate and textThemes and motifsCriticism and interpretationLegacyScene by sceneSee alsoNotes and referencesSourcesExternal linksNavigation menu"Consumer Price Index (estimate) 1800–"10.2307/28710160037-3222287101610.1093/res/II.5.31910.2307/45967845967810.2307/2869925286992510.1525/jams.1982.35.3.03a00050"Dada Masilo: South African dancer who breaks the rules"10.1093/res/os-XV.57.1610.2307/28680942868094"Sweet Sorrow: Mann-Korman's Romeo and Juliet Closes Sept. 5 at MN's Ordway"the original10.2307/45957745957710.1017/CCOL0521570476.009"Ram Leela box office collections hit massive Rs 100 crore, pulverises prediction"Archived"Broadway Revival of Romeo and Juliet, Starring Orlando Bloom and Condola Rashad, Will Close Dec. 8"Archived10.1075/jhp.7.1.04hon"Wherefore art thou, Romeo? To make us laugh at Navy Pier"the original10.1093/gmo/9781561592630.article.O006772"Ram-leela Review Roundup: Critics Hail Film as Best Adaptation of Romeo and Juliet"Archived10.2307/31946310047-77293194631"Romeo and Juliet get Twitter treatment""Juliet's Nurse by Lois Leveen""Romeo and Juliet: Orlando Bloom's Broadway Debut Released in Theaters for Valentine's Day"Archived"Romeo and Juliet Has No Balcony"10.1093/gmo/9781561592630.article.O00778110.2307/2867423286742310.1076/enst.82.2.115.959510.1080/00138380601042675"A plague o' both your houses: error in GCSE exam paper forces apology""Juliet of the Five O'Clock Shadow, and Other Wonders"10.2307/33912430027-4321339124310.2307/28487440038-7134284874410.2307/29123140149-661129123144728341M"Weekender Guide: Shakespeare on The Drive""balcony"UK public library membership"romeo"UK public library membership10.1017/CCOL9780521844291"Post-Zionist Critique on Israel and the Palestinians Part III: Popular Culture"10.2307/25379071533-86140377-919X2537907"Capulets and Montagues: UK exam board admit mixing names up in Romeo and Juliet paper"Istoria Novellamente Ritrovata di Due Nobili Amanti2027/mdp.390150822329610820-750X"GCSE exam error: Board accidentally rewrites Shakespeare"10.2307/29176390149-66112917639"Exam board apologises after error in English GCSE paper which confused characters in Shakespeare's Romeo and Juliet""From Mariotto and Ganozza to Romeo and Guilietta: Metamorphoses of a Renaissance Tale"10.2307/37323537323510.2307/2867455286745510.2307/28678912867891"10 Questions for Taylor Swift"10.2307/28680922868092"Haymarket Theatre""The Zeffirelli Way: Revealing Talk by Florentine Director""Michael Smuin: 1938-2007 / Prolific dance director had showy career"The Life and Art of Edwin BoothRomeo and JulietRomeo and JulietRomeo and JulietRomeo and JulietEasy Read Romeo and JulietRomeo and Julieteeecb12003684p(data)4099369-3n8211610759dbe00d-a9e2-41a3-b2c1-977dd692899302814385X313670221313670221