I have a puzzling problem using GeoTiff rasters which have been generated using Python + GDAL functions. When these grids are "Sampled" (ArcMap 10.4.1 Spatial Analyst Sample command), some rows of data contain either Nodata values, or in some cases very small or large negative values where there should be valid data values.
This occurs when sampling a grid, using the grid itself as the location of points to be sampled. There are no errors or warnings reported by Sample.
This image shows the X, Y values output for a grid (using Sample tool) which has been loaded into ArcMap with the values symbolised to show NoData (-9999) values in red. The output points all have the correct coordinates, and there are not gaps apparent. However, when symbolised to the value of the underlying raster, you can see the values in two lines across the grid either appear with NoData values, or in the case of the lower line, contain values which are not the same as the surrounding grid values that would be expected. Inspecting the individual values in those rows, none of the extracted values correspond to any of the grids being sampled at that point. The lines are 140 m apart (70 rows of 2m grid)
Initially I was sampling several grids at the same time, but I found that I could reproduce the problem using only one of those grids at a time. I have tried many things to find a cause/solution this including:
- defining the projection (the same as the grid already has defined),
- sample using a similar grid created in ArcMap with same grid size,
- exporting a world file for the Geotiff.
but these do not work.
Two things that DO have an effect are
- converting the grid to a point Shapefile, and using the point Shapefile to sample with, instead of the raster.
- I also found that projecting the raster (from GDA_1994_MGA_Zone_55 to GDA_1994_MGA_Zone_55) changes the spacing of the invalid rows. In the case of the grid shown, the rows disappear, but with a larger grid, the spacing of invalid rows simply increased
This suggests there is something going on in ArcMap's interpretation of the grid projection, which was based on EPSG in the GDAL functions. I am wondering if anyone has seen this issue with ArcGIS 10.x, and would like to know if there is a way I can resolve this in the Python GDAL functions used to generate the grid (or perhaps is this an ArcMap bug?).
