Wednesday, August 31, 2011

Latitude & Longitude rasters in ArcGIS 10

The ability to calculate latitude and longitude rasters using the GRID commands $$YMap and $$XMap was removed in ArcGIS 10 (as well as all the GRID commands including $$COLMAP, $$ROWMAP, $$NROWS, $$NCOLS, $$CELLSIZE). To get around this you can still load the ArcGIS 9.3 geoprocessor in ArcGIS 10 and then call the functions using the Single Output Map Algebra tool. Latitude & longitude rasters can also be calculated with NumPy or the FlowAccumulation function but these aren't great options for really large rasters. Make sure to use a DEM that is in a geographic coordinates system and then project the results.
def lat_long_rasters(ned_raster, lat_raster, lon_raster):
    import arcgisscripting
    gp = arcgisscripting.create(9.3)
    gp.CheckOutExtension("Spatial")
    gp.WorkSpace = os.path.dirname(ned_raster)
    gp.SingleOutputMapAlgebra_sa(
        "Con(IsNull({0}), {1}, $$YMap)".format(ned_raster, ned_raster), 
        lat_raster)
    gp.SingleOutputMapAlgebra_sa(
        "Con(IsNull({0}), {1}, $$XMap)".format(ned_raster, ned_raster), 
        lon_raster)
    del gp

Friday, August 19, 2011

Nan in NumPy array to NoData in ArcGIS raster

For smaller rasters, NumPy is a much faster alternative to some of the ArcGIS geoprocessing functions. One issue with doing this is that every raster must have exactly the same dimensions (rows, columns, cellsize, extent, etc.) since NumPy has no spatial component. Later on I will post about some of the techniques we use to make sure they are all identical.

There is one issue with converting NumPy arrays back to rasters though that I want to outline. The arcpy function NumPyArrayToRaster does not correctly convert NaN values. The values will display as NoData in the identity window, but they do not act as NoData in subsequent calculations and the 'Display NoData as' option in the Symbology tab has no impact.

To get around this, convert all NaN values to some unused value, and then set that value as the value_to_nodata parameter. For example:

test_array[np.isnan(test_array)] = -9999
test_obj = arcpy.NumPyArrayToRaster(test_array, pnt, cellsize, cellsize, -9999)
test_obj.save(test_raster)

This is Bug NIM063028 and is not fixed in any current or future version.

Pyramids documentation

There is a small issue with the documentation for the pyramids environment parameter. The explanation section for the pyramid_option has a few errors. The option is "PYRAMIDS", not "PYRAMID", and the default is to build full pyramids (-1 level), not NONE. This is Bug NIM062605 and is reported as fixed in 10.1. The example code at the bottom of the help page shows the correct usage. To force no pyramids I prefer to use: env.pyramid = "PYRAMIDS 0".

Also, the CreateRandomRaster function doesn't honor the pyramids setting. This is Bug NIM062606 and is not currently reported as fixed.

Wednesday, August 17, 2011

Offset of TIF rasters in ArcGIS 10, pre-SP2

In Arc10, pre-SP2, our Landsat TIF rasters were being shifted half a pixel away from where they should be.  For older Landsat 5TM imagery the thermal band could be up 60m offset in both the x and y.  ESRI fixed the problem in SP2 (Bug NIM060759) but in case you are loading TIF rasters with a pre-SP2 version of Arc10, it is easy to shift the rasters back to the correct location.  Since we were not using the panchromatic rasters from Landsat 7, all of our half cell sizes were integers (30m -> 15m, 60m -> 30m, 120m -> 60m), but that may be an issue for other applications.
half_cell_width = int(0.5*arcpy.Describe(landsat_tif).meanCellWidth+0.5)
half_cell_height = int(0.5*arcpy.Describe(landsat_tif).meanCellHeight+0.5)
arcpy.Shift_management(landsat_tif, shift_raster, -half_cell_width, half_cell_height)


It is also fairly simple to check which ArcGIS version and service pack you have:
version_str = arcpy.GetInstallInfo("desktop")["SPBuild"]
version_str = "10.0.2.3200"
sp_str = version_str.split(".")[2]
sp_str = "2"

MakeFeatureLayer lock files

When calling the MakeFeatureLayer function in a python script, a lock file is created that may persist even after the script finishes.  To remove the lock file, just call arcpy.Delete_management(layer) before the script terminates.

This forum post provided the answer to this problem.