In this blog post, I want to get the Ever Given cargo ship stuck in the pond next to MathWorks at scale. This blog post is inspired by Ned Gulley using Ever Given Ever Ywhere.
We'll download the ship from a public source and isolate it.
ship = imread('https://e3.365dm.com/21/03/768x432/skynews-ever-given-suez-canal_5320248.jpg?20210327160438'); imshow(ship)
I could probably come up with a way to automatically segment the ship. But I only need to do it once, so I'll just draw the polygon and mask it.
if false % Draw polygon on border of ship p = drawpolygon; else load('p.mat', "p") imshow("eg.png") end
Mask the ship based on the polygon and apply that to the original image.
shipmask = poly2mask(p.Position(:, 1), p.Position(:, 2), height(ship), width(ship)); shipmask = imerode(shipmask, strel("disk", 2)); figure tiledlayout(2, 1, TileSpacing="compact", Padding="tight") nexttile imshow(shipmask) nexttile ship = ship.*cast(shipmask, class(ship)); imshow(ship)
Get the orientation from the ship mask and apply it to the original image.
shipprops = regionprops(shipmask, "Orientation") shiphorizontal = imrotate(ship, -shipprops.Orientation); figure imshow(shiphorizontal)
shipprops = struct with fields: Orientation: 19.588174604138796
Crop to just the ship.
shipbbox = regionprops(logical(shiphorizontal(:, :, 1)), "BoundingBox") shipthumb = imcrop(shiphorizontal, shipbbox.BoundingBox); imshow(shipthumb)
shipbbox = struct with fields: BoundingBox: [1.815000000000000e+02 3.205000000000000e+02 427 67]
We'll download the satellite imagery from USGS for the region we care about. I'll interactively pick the regional limits for the area I care about.
if false % Interactively select limits w = webmap [latlim, lonlim] = wmlimits(w) else latlim = [42.2966508472710 42.3053796253471] lonlim = [-71.3800315706306 -71.3639383165411] end
latlim = 42.296650847271003 42.305379625347101 lonlim = -71.380031570630607 -71.363938316541095
Find the ortho layer and download it for this region.
ortho = wmsfind('usgsimageryonly', SearchField='serverurl'); [lakeside, lakesidereference] = wmsread(ortho, latlim=latlim, lonlim=lonlim); figure geoshow(lakeside, lakesidereference) axis tight
I want to put the ship in the little lake between Rt 9 and the railroad track immediately adjacent to MathWorks' Lakeside campus. It looks like a ship orientation of about will work well.
shipn75 = imrotate(shipthumb, -75,'nearest'); figure imshow(shipn75);
We'll work in meters and need the dimensions of the ship. A quick web search says the length is 399.94m and the beam is 58.8m.
shiplength = 399.94; shipbeam = 58.8;
At the same time will check the ratio of length to width with that of our image.
props = regionprops(logical(shipn75(:,:,1)), "MajorAxisLength", "MinorAxisLength"); image_shipratio = props.MajorAxisLength/props.MinorAxisLength actual_shipratio = shiplength/shipbeam
image_shipratio = 6.743610859446426 actual_shipratio = 6.801700680272109
I'm happy with that(!) especially considering I just approximated the ship with a drawn polygon.
Now we need to figure out the size of this rotated image in meters. Assuming the ship was a perfect rectangle, we can calculate this size using some trigonometry. It's been awhile since I've done this so I broke out a trusty envelope and tried to remember soh-cah-toa. The equation for height is:
imgheightm = shiplength+sind(15)*shipbeam
imgheightm = 4.151585598520282e+02
I can use the ratio here to figure out the width in meters since the pixels are isotropic.
imgwidthm = imgheightm/height(shipn75)*width(shipn75)
imgwidthm = 1.700996877171504e+02
Finally, we need to convert the background of the image to NaN so when we plot it shows as transparent. This is where all channels are 0.
shipoverlay = im2double(shipn75); shipoverlay(repmat(all(~shipoverlay, 3), [1 1 3])) = NaN;
We need to project the satellite imagery from geographic coordinates (lat/lon) to cartesian coordinates (x/y in meters) to overlay the boat on it. We'll use the Mass State Plane coordinate reference system which is optimized for Massachusetts. A quick web search tells me that the code for this coordinate system is 26986.
p = projcrs(26986) [latgrid, longrid] = geographicGrid(lakesidereference); [xgrid, ygrid] = projfwd(p, latgrid, longrid); figure mapshow(xgrid,ygrid,lakeside) axis tight
p = projcrs with properties: Name: "NAD83 / Massachusetts Mainland" GeographicCRS: [1×1 geocrs] ProjectionMethod: "Lambert Conic Conformal (2SP)" LengthUnit: "meter" ProjectionParameters: [1×1 map.crs.ProjectionParameters]
I'll use a datatip to find the upper left corner of where the boat should go. The live editor gives me the below code to
ax = gca; chart = ax.Children(1); datatip(chart, 210474, 894667);
Since the map is now in projected coordinates, the x/y values represent meters. We can create map reference cells that describe the pixel size of the boat image and location which then makes plotting it easy.
xpt = 210474; ypt = 894667; mr = maprefcells([xpt xpt+imgwidthm], [ypt-imgheightm ypt], size(shipoverlay, [1 2]), ColumnsStartFrom="north")
mr = MapCellsReference with properties: XWorldLimits: [210474 210644.099687717] YWorldLimits: [894251.841440148 894667] RasterSize: [432 177] RasterInterpretation: 'cells' ColumnsStartFrom: 'north' RowsStartFrom: 'west' CellExtentInWorldX: 0.961015184842692 CellExtentInWorldY: 0.961015184842701 RasterExtentInWorldX: 170.099687717156 RasterExtentInWorldY: 415.158559852047 XIntrinsicLimits: [0.5 177.5] YIntrinsicLimits: [0.5 432.5] TransformationType: 'rectilinear' CoordinateSystemType: 'planar' ProjectedCRS: 
And now we can finally plot the Lakeside image with the Ever Given overlaid on top to scale.
figure hold on mapshow(xgrid, ygrid, lakeside) s = mapshow(zeros(size(shipoverlay, [1 2])), mr, DisplayType="surface", CData=shipoverlay); axis tight xticks() yticks()
Published with MATLAB® R2021a
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.