search

zoran-cuckovic / QGIS-visibility-analysis

Quantum GIS plugin for visibility analysis
GNU General Public License v3.0
73 stars 18 forks source link

Geolocation does not seem correct in output #4

Closed xrubio closed 7 years ago

xrubio commented 7 years ago

I have been trying to use the plugin to generate a viewshed for a point layer in Britain but it seems that the result is not correct (see attached options and result in the screenshots). It can be seen how the densest zones are several kilometers apart from where they should be. The CRS is epsg:27700 (OSGB 1936 / British National Grid).

I am using QGIS 2.16.3 but I also tested it with 2.18 and 2.8 Wien and I get the same results.

options result

zoran-cuckovic commented 7 years ago

Hello,

This has to be an issue of coordinate reference system. Try to switch off the automatic "on the fly" re-projection and see whether everything remains on its proper position (points, DEM and the result). The plugin does not make any adjustment of coordinate systems on its own.

zoran-cuckovic commented 7 years ago

PS a wild guess : your DEM is in WGS 84 (judging by elliptic viewsheds)?

xrubio commented 7 years ago

Hi Zoran and thank you for creating the plugin!

Everything is correct in the CRS of my project (all layers have the same CRS so on-the-fly do not change anything). In addition I also tried the plugin with a different project in Greece (epsg 32634, WGS 84 UTM 34N) and I got the same issue.

The ellipsoid is airy:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

I tried changing the 'ellipsoid for distance calculations' to no effect. I would say that the x coordinates are correct and there is something going on with the Y axis (I attach image).

Xavi

shift

xrubio commented 7 years ago

Quick update...the plugin works perfect with the Greece project. Could it be something related to the British CRS or the airis ellipsoid...?

zoran-cuckovic commented 7 years ago

Well, if anything, elliptic viewsheds cannot happen unless due to a re-projection of the output raster. I cannot imagine another type of issue.

Are you willing to share a sample of your data (just a few points and a part of the DEM)? I could then try for myself.

xrubio commented 7 years ago

Sure! I uploaded a zipe file with the clipped DEM and the shapefile here: https://we.tl/dqcMWKLoe2

2016-11-07 18:35 GMT+00:00 Zoran Čučković notifications@github.com:

Well, if anything, elliptic viewsheds cannot happen unless due to a re-projection of the output raster. I cannot imagine another type of issue.

Are you willing to share a sample of your data (just a few points and a part of the DEM)? I could then try for myself.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/zoran-cuckovic/QGIS-visibility-analysis/issues/4#issuecomment-258921545, or mute the thread https://github.com/notifications/unsubscribe-auth/ABWcSrTPh9hBlglDzcCvHurBc33IlN4Lks5q729bgaJpZM4KrBYh .

zoran-cuckovic commented 7 years ago

Here is an example of elliptic output : the data is SRTM elevation model, in WGS, but reprojected "on the fly" to your CRS (OSGB 1936 / British National Grid).

scotland_viewsheds

And considering the drift of several hundred kilometres to the south: I would say your points have a different CRS than the DEM.

You need to re-project your DEM into an appropriate CRS prior to analysis. Then you need to make sure your points are in the same CRS.

xrubio commented 7 years ago

But the shapefiles and the DEM are already reprojected to OSGB (I used gdal_warp with DEMs from SRTM). You can see here the output of gdalinfo:

gdalinfo dem.tif

Driver: GTiff/GeoTIFF

Files: dem.tif

dem.tif.aux.xml

Size is 19938, 12349

Coordinate System is:

PROJCS["OSGB 1936 / British National Grid",

GEOGCS["OSGB

1936",

DATUM["OSGB_1936",

        SPHEROID["Airy

1830",6377563.396,299.3249646,

AUTHORITY["EPSG","7001"]],

TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],

AUTHORITY["EPSG","6277"]],

PRIMEM["Greenwich",0,

        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4277"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",49],
PARAMETER["central_meridian",-2],
PARAMETER["scale_factor",0.9996012717],
PARAMETER["false_easting",400000],
PARAMETER["false_northing",-100000],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","27700"]]

Origin = (-172964.868743266095407,1150958.007875331910327) Pixel Size = (53.892136284530039,-92.307591020375554) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( -172964.869, 1150958.008) ( 12d15' 8.76"W, 59d50'37.71"N) Lower Left ( -172964.869, 11051.566) ( 9d57'16.56"W, 49d43'30.75"N) Upper Right ( 901536.544, 1150958.008) ( 6d59'38.06"E, 59d56'11.66"N) Lower Right ( 901536.544, 11051.566) ( 4d58'16.44"E, 49d47'20.08"N) Center ( 364285.838, 581004.787) ( 2d33'36.24"W, 55d 7'19.49"N) Band 1 Block=19938x1 Type=Int16, ColorInterp=Gray Min=-17.000 Max=1211.000 Minimum=-17.000, Maximum=1211.000, Mean=137.169, StdDev=142.255 NoData Value=-32768 Metadata: STATISTICS_MAXIMUM=1211 STATISTICS_MEAN=137.16900903408 STATISTICS_MINIMUM=-17 STATISTICS_STDDEV=142.25489207991

2016-11-07 19:08 GMT+00:00 Zoran Čučković notifications@github.com:

Here is an example of elliptic output : the data is SRTM elevation model, in WGS, but reprojected "on the fly" to your CRS (OSGB 1936 / British National Grid).

[image: scotland_viewsheds] https://cloud.githubusercontent.com/assets/6622934/20071281/20715be8-a524-11e6-89eb-4c1a608a1254.jpg

And considering the drift of several hundred kilometres to the south: I would say your points have a different CRS than the DEM.

You need to re-project your DEM into an appropriate CRS prior to analysis. Then you need to make sure your points are in the same CRS.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/zoran-cuckovic/QGIS-visibility-analysis/issues/4#issuecomment-258931395, or mute the thread https://github.com/notifications/unsubscribe-auth/ABWcSpNH7G6uw8Xmua5EoPN1G9dqQnAOks5q73dKgaJpZM4KrBYh .

zoran-cuckovic commented 7 years ago

The problem is that you have rectangular pixels (!) Pixel Size: 53.8914,-92.3058

pixels I've never had to do with such a dataset. You need to reproject to a regular pixel grid. See eg. : http://gis.stackexchange.com/questions/27924/how-to-square-the-pixels-of-a-geotiff

xrubio commented 7 years ago

Yeah I suspected that this could be the issue after checking gdalinfo....ok, I'll try this solution, thanks for your help!

Xavi

2016-11-07 19:24 GMT+00:00 Zoran Čučković notifications@github.com:

The problem is that you have rectangular pixels (!) Pixel Size: 53.8914,-92.3058

[image: pixels] https://cloud.githubusercontent.com/assets/6622934/20072175/aec9212a-a527-11e6-8bcf-2e24a2bf0dd0.jpg I've never had to do with such a dataset. You need to reproject to a regular pixel grid. See eg. : http://gis.stackexchange.com/ questions/27924/how-to-square-the-pixels-of-a-geotiff

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/zoran-cuckovic/QGIS-visibility-analysis/issues/4#issuecomment-258935779, or mute the thread https://github.com/notifications/unsubscribe-auth/ABWcSi9H2wWJFTGPTtF-Jp_TdEJBThVIks5q73sBgaJpZM4KrBYh .

xrubio commented 7 years ago

Ok it worked like a charm, so the problem was that pixels were non-squared.