How to plot a 3d surface with a geotiff on top?

See below for my code.

I have a side scan sonar that kind of sucks, the data is proprietary but I wrote a python script that will read and generate a geotiff of the sonar. I have a matlab script that will read the single beam data and generate a 3d surface of the bathymetry, what I want it to be able to plot the geotiff Im creating on top of the 3d surface so you can see the sidescan but also see depth information. I can do both of them seperatly but I cant get them together. Any ideas on getting a geotiff on top of a 3d surface? Thanks!

P1 = uigetdir('C:\');
S10 = dir(fullfile(P1,"*.tif"));
F11 = fullfile(S10.folder,S10.name);
[A, R, cmap] = readgeoraster(F11);

depth = -inst_dep_m;

LatLong = [lat, lon, depth];
[MM, ~, ~] = rmoutliers(LatLong, 1);
lat1 = MM(:,1);
lon1 = MM(:,2);
depth1 = MM(:,3);

depth1(end) = NaN;
a = depth1;

[xData, yData, zData] = prepareSurfaceData( lon1, lat1, depth1 );

ft = fittype( 'lowess' );
[fitresult, ~] = fit( [xData, yData], zData, ft, 'Normalize', 'on' );

figure('Name','True Bathy Surface', 'Position',[850, 125, 600, 550]);
ax1 = axes();
mapshow(A, cmap, R)
plot(fitresult);
hold on 
xlabel('Longitude (°)', 'Interpreter', 'none' );
ylabel('Latitude (°)', 'Interpreter', 'none' );
zlabel('Depth (ft)', 'Interpreter', 'none' );
hold on 
h1 = colorbar('southoutside');
h1.Label.String = 'Depth (ft)';
colormap
view(3)
view(30,65);
rotate3d on;

% colormap(ax1, 'turbo');

[Walter Roberson]

The "hold on" might need to be before the plot(fitresult)

You should not need two "hold on"

[Bradley]

Im not sure why I had so many hold on in there. This is what I changed it too and nothing happened

plot(fitresult);
hold on 
mapshow(A, cmap, R)
xlabel('Longitude (°)', 'Interpreter', 'none' );
ylabel('Latitude (°)', 'Interpreter', 'none' );
zlabel('Depth (ft)', 'Interpreter', 'none' );

I cant tell what is supposed to be happening. I get the same surface I had before which is good that its still being formed but it doesnt include the geotiff on that surface, which is what I want. Thanks!

Answer 1

[Umar]

Hi @Bradley ,

To achieve the desired overlay of your geotiff on the 3D surface plot, you need to ensure that both datasets are aligned correctly in terms of their coordinate systems. Here’s a structured approach to help you integrate these visualizations effectively:

1. Coordinate System Alignment: Ensure that the coordinates used in your geotiff match those of your bathymetric surface. If the geotiff is not aligned with the bathymetric data, it will not display correctly on top of it.

2. Modify Your Code: You can adjust your plotting sequence and ensure you're using the correct axes for each plot. Here's an updated version of your code snippet that should help:

     % Load Geotiff
      P1 = uigetdir('C:\');
      S10 = dir(fullfile(P1,"*.tif"));
      F11 = fullfile(S10.folder,S10.name);
      [A, R, cmap] = readgeoraster(F11);

      % Prepare Depth Data
      depth = -inst_dep_m;
      LatLong = [lat, lon, depth];
      [MM, ~, ~] = rmoutliers(LatLong, 1);
      lat1 = MM(:,1);
      lon1 = MM(:,2);
      depth1 = MM(:,3);
      depth1(end) = NaN;

      % Prepare Surface Data
      [xData, yData, zData] = prepareSurfaceData(lon1, lat1, depth1);
      ft = fittype('lowess');
      [fitresult, ~] = fit([xData, yData], zData, ft, 'Normalize', 'on');

      % Create Figure
      figure('Name','True Bathy Surface', 'Position',[850, 125, 600, 550]);
      ax1 = axes();

      % Plotting 3D Surface
      plot(fitresult);
      hold on;

      % Overlay Geotiff
      mapshow(A, cmap, R); 
     % Make sure R corresponds to the correct geographic reference
     xlabel('Longitude (°)', 'Interpreter', 'none');
     ylabel('Latitude (°)', 'Interpreter', 'none');
     zlabel('Depth (ft)', 'Interpreter', 'none');

      % Colorbar and View Settings
      h1 = colorbar('southoutside');
      h1.Label.String = 'Depth (ft)';
      colormap(ax1, 'turbo');

      view(3);
      rotate3d on;

      % Adjust transparency if needed
      alpha(0.5); % Adjust transparency to see through the geotiff if required

3. Key Adjustments:

Order of Plotting: The `mapshow` command should be placed after plotting your fit result but before setting labels or modifying views.

Transparency: If you want to visualize both layers clearly without one obscuring the other, consider adjusting the alpha transparency of your geotiff overlay.

Axes Limits: Ensure that your axes limits are set appropriately to encompass both datasets. This can be done using `xlim`, `ylim`, and `zlim` commands.

Georeferencing: If your geotiff is not displaying correctly, double-check its georeferencing information (the variable `R`). The geographic referencing must match that of your bathymetric data for accurate overlay.

Debugging: If issues persist after these adjustments, consider isolating each plotting step to verify that each layer appears correctly before integrating them.

Visualization Techniques: Depending on the complexity of your datasets and intended analysis, you might explore additional visualization techniques or libraries such as `surf` or `meshgrid` for more advanced surface rendering.

By following these guidelines and refining your code accordingly, you should be able to successfully overlay your geotiff onto the 3D surface plot in MATLAB.

Hope this helps.