Problem with using viewshed in different geographical locations.

viewshed produces invalid visibility matrix that is all zeroes except for one entry when given SRTM30+ geographical raster data that is near latitude 25 deg N, longitude 54 deg E. Running the same script with data near latitude 23 deg N, longitude 121 deg E works fine. GeographicCellsReference with properties data for the case that works is

GeographicCellsReference with properties:

LatitudeLimits: [10.0041666666667 45.0041666666667]

LongitudeLimits: [100.004166666667 140.004166666667]

RasterSize: [4201 4801]

RasterInterpretation: 'cells'

ColumnsStartFrom: 'north'

RowsStartFrom: 'west'

CellExtentInLatitude: 0.00833134967864794

CellExtentInLongitude: 0.0083315975838367

RasterExtentInLatitude: 35

RasterExtentInLongitude: 40

XIntrinsicLimits: [0.5 4801.5]

YIntrinsicLimits: [0.5 4201.5]

CoordinateSystemType: 'geographic'

GeographicCRS: []

AngleUnit: 'degree'

GeographicCellsReference with properties data for the case that does not work is

GeographicCellsReference with properties:

LatitudeLimits: [10.0041666666667 38.0041666666667]

LongitudeLimits: [40.0041666666667 74.0041666666667]

RasterSize: [3361 4081]

RasterInterpretation: 'cells'

ColumnsStartFrom: 'north'

RowsStartFrom: 'west'

CellExtentInLatitude: 0.00833085391252603

CellExtentInLongitude: 0.00833129135015927

RasterExtentInLatitude: 28

RasterExtentInLongitude: 34

XIntrinsicLimits: [0.5 4081.5]

YIntrinsicLimits: [0.5 3361.5]

CoordinateSystemType: 'geographic'

GeographicCRS: []

AngleUnit: 'degree'

MSL elevation data provided to viewshed was all double precision, real, positive, greater than or equal to zero with no nan entries for both cases.Transmit location was within the coverage region for both cases. Command used for visibility was:

vis = viewshed(Zterrain,Rterrain,rdrlat,rdrlon,rdralt,tgtalt,"MSL","MSL",Re,5/2*Re);

with

Re = 6371000

Zterrain = terrain raster 2D array (meters above MSL)

Rterrain = GeographicCellsReference with properties

rdralt = observer location altitude scalar variable (meters)

Plots using commands:

figure

mapshow(Zterrain,Rterrain,'DisplayType','mesh')

demcmap(Zterrain)

produce valid results.

Answer 1

[Divit]

Hi Paul,

Given the information you've provided, here are a few aspects to consider that might help resolve the issue:

• Observer and Target Altitudes: Ensure that the observer's altitude ("rdralt") and the target altitude ("tgtalt") are correctly specified and that they are within reasonable ranges for the given terrain data. If the observer's altitude is not significantly higher than the surrounding terrain, the viewshed could be obscured, leading to an almost entirely zero visibility matrix.
• Observer Location: Verify that the observer location ("rdrlat", "rdrlon") is within the bounds of the raster data and not on the edge, which could cause issues with visibility calculations.
• Earth's Radius (Re): Using a constant value for the Earth's radius is generally okay, but if you need more precision, consider using a more accurate model that accounts for the Earth's ellipsoidal shape.
• Coordinate Reference System (CRS): The "GeographicCRS" property is empty in both cases. If your data requires a specific CRS, make sure to define it.

If after checking all these factors the problem persists, consider the following troubleshooting steps:

• Subsample the Data: Try using a subset of your raster data that is closer to the observer location to see if the issue is related to the extent of the data.
• Adjust Observer Parameters: Temporarily modify the observer parameters (altitude, location) to different values to see if there's a threshold where the viewshed starts producing non-zero results.

You can read more about "viewshed" by referrring to the following documentation link:

https://www.mathworks.com/help/map/ref/viewshed.html