2

Is there a possibility to calculate the horizon line of a point with global mapper, not the viewshed, but just a line (or polygon) joining all the furthermost points from A and are at the same time seen from it.

  • I found the scripts I used, so I just updated my answer, adding the scripts to it. I also added a few figures of a example I just run. – Camilo Rada Sep 11 '18 at 1:06
2

No, it is not possible. However I've done it with a workaround that is a bit convoluted but works.

Path profile is one of the Global Mapper features that can be scripted in Global Mapper. So you can create an script that calculates path profiles from from your point of view all around at your desired resolution (I did it every 0.5° I think), and exports the results to .xyz files (horizontal coordinates and elevation of each point).

Then I created an script (in Matlab, but Python or anything similar would work), that reads all the .xyz files and finds the point at higher elevation (angle above the horizon from the point of view) in each file and write it in a new output file. Therefore, that file contains exactly what you are looking for. You can easily import that file back into Global Mapper and plot it as a horizon line.

Scripting makes sense if you will do it many times. For just one time it might be quicker to just build the line manually based on the view shed output.

The first block of the Global Mapper script would look like this

GLOBAL_MAPPER_SCRIPT VERSION="1.00"
GENERATE_PATH_PROFILE \
    FILENAME="C:\LOS\LOS_000_deg.xyz" \
    ELEV_UNITS=METERS \
    POINT_COUNT=1000 \
    START_POS=334096.311,6316949.632 \
    END_POS=293308.858,6343374.188 \

And you can easily make another script to generate an script like the one above but with one GENERATE_PATH_PROFILE block for each azimuth you want to find an horizon point for.

If you are familiar with Matlab you can use my scripts. The first one is

E0=334000;
N0=6317000;
Dist=100000;
dist_resolution=15;
azimuth_resolution=0.5;
outfile='GM_script.gmw';

fid=fopen(outfile,'w');
fprintf(fid,'GLOBAL_MAPPER_SCRIPT VERSION="1.00"\n');

for theta=0:azimuth_resolution:360-azimuth_resolution
    DN=sind(theta)*Dist;
    N=N0+DN;
    DE=cosd(theta)*Dist;
    E=E0+DE;
    fprintf(fid,'GENERATE_PATH_PROFILE FILENAME="C:\\users\\camilo\\My Documents\\Desktop\\tmp\\LOS\\LOS_%06d_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=%d START_POS=%f,%f END_POS=%f,%f\n',theta*1000,round(Dist/dist_resolution),E0,N0,E,N);
end

fclose(fid);

Where E0 and N0 are the easting and northing of the point you are calculating the horizon from. dist_resolution is resolution in meters along each path profile. azimuth_resolution the one across the horizon in degrees. And Dist the maximum distance to search horizon points.

This generates the script GM_script.gmw, that looks like this:

GLOBAL_MAPPER_SCRIPT VERSION="1.00"
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_000000_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=434000.000000,6317000.000000
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_000500_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=433996.192306,6317872.653550
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_001000_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=433984.769516,6318745.240644
...
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_358500_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=433965.732498,6314382.305169
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_359000_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=433984.769516,6315254.759356
GENERATE_PATH_PROFILE FILENAME="C:\users\camilo\My Documents\Desktop\tmp\LOS\LOS_359500_mdeg.xyz" ELEV_UNITS=METERS POINT_COUNT=6667 START_POS=334000.000000,6317000.000000 END_POS=433996.192306,6316127.346450

Then you have to run GM_script.gmw in Global Mapper from the File menu "Run Script..." like this (make sure to check the context box)

enter image description here

In this example, that will generate 720 .xyz files. The script that compile the data of all files is as follow

outfile='horizon_1.csv';
E0=334000;
N0=6317000;
H0=493.659+1.5;

files=dir('*_mdeg.xyz');
nfiles=length(files);
fid=fopen(outfile,'w');
horizon=zeros(nfiles,5);
for i=1:nfiles
    disp(['Reading ' files(i).name '...'])
    data=load(files(i).name);
    E=data(2:end,1);
    N=data(2:end,2);
    H=data(2:end,3);
    D=sqrt(((E-E0).^2)+((N-N0).^2));
    elev=atand((H-H0)./D);
    [~,idx]=max(elev);
    fprintf(fid,'%.3f,%.3f,%.3f\n',E(idx),N(idx),H(idx));
    theta=str2double(files(i).name(5:10))/1000;
    horizon(i,:)=[E(idx),N(idx),H(idx),elev(idx),theta];
end
fclose(fid);
figure
plot(horizon(:,5),horizon(:,4));
xlim([0,360])
xlabel('Azimuth [°]')
ylabel('Horizon elevation [°]')

This creates the horizon_1.csv file that if loaded in Global Mapper you will get the line you are looking for. In this case it looks like this

enter image description here

The script also generates a plot with the profile of the horizon. Like this one

enter image description here

With an slightly more sophisticated script you can also plot all the visible points color-coded according to their distance. Therefore you can see how the different areas of the terrain overlay each other and how far away is the horizon on each direction. Like this

enter image description here

  • Thank you @Camilo so much, I appreciate your work, that is so useful, I am familiar wìth QGIS Desktop, I think I would ask the same question but with Qgis, thank you again – Khaled Sep 11 '18 at 10:39
  • @ camilo could you please send me your "GM_script.gmw" file? – Khaled Sep 13 '18 at 10:09
  • @Khaled Sure, here is the one I used in the example above: dropbox.com/s/qevnbbyl14al3w9/GM_script.gmw?dl=0 note that commands in Global Mapper scripts go in one line, but you can split them in multiple lines using "\" as I did in the first example above (for readability). If you are not familiar with matlab scripts you can easily generate a similar script file using Excel or OpenOffice Calc or any spreadsheet. I also added and example in my answer. – Camilo Rada Sep 13 '18 at 16:40
  • Believe me, I don't know how to think you, if you allow, I would ask also sthg :( , if I use the geographical coordinates lon/lat for my points instead of Easting/Northing, could the $ELEV_UNITS=METERS$ and $POINT_SPACING=10.0$ remain working? bcz I don't know so much about UTM specially the area I am interesting for, is a small country 2 degree wide and 2 degree height on map, but there is 2 UTM zones overlapping :( – Khaled Sep 13 '18 at 17:01
  • @Khaled I think you can used whatever coordinates are in the workspace you are running the script over. So if the workspace is lat/lon, you can use lat/lon. Elev_units and spacing should remain the same. – Camilo Rada Sep 13 '18 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.