{"id":12187,"date":"2021-04-02T16:08:08","date_gmt":"2021-04-02T20:08:08","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=12187"},"modified":"2021-04-10T17:30:44","modified_gmt":"2021-04-10T21:30:44","slug":"get-the-ship-stuck","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2021\/04\/02\/get-the-ship-stuck\/","title":{"rendered":"Get the Ship Stuck!"},"content":{"rendered":"<div class=\"content\">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 <a href=\"https:\/\/evergiven-everywhere.glitch.me\/\">Ever Given Ever Ywhere<\/a>.&nbsp;<\/p>\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li><a href=\"#1\">Get the Ship<\/a><\/li>\n<li><a href=\"#6\">Download Lakeside Satellite Imagery<\/a><\/li>\n<li><a href=\"#8\">Prepare the Ship for \"Docking\"<\/a><\/li>\n<li><a href=\"#13\">Prepare Lakeside for a Ship<\/a><\/li>\n<li><a href=\"#16\">Stick the Ship!<\/a><\/li>\n<\/ul>\n<\/div>\n<h3>Get the Ship<a name=\"1\"><\/a><\/h3>\n<p>We'll download the ship from a public source and isolate it.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">ship = imread(<span style=\"color: #a020f0;\">'https:\/\/e3.365dm.com\/21\/03\/768x432\/skynews-ever-given-suez-canal_5320248.jpg?20210327160438'<\/span>);\r\nimshow(ship)<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_01.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>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.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\"><span style=\"color: #0000ff;\">if<\/span> false\r\n<span style=\"color: #228b22;\">% Draw polygon on border of ship<\/span>\r\n    p = drawpolygon;\r\n<span style=\"color: #0000ff;\">else<\/span>\r\n    load(<span style=\"color: #a020f0;\">'p.mat'<\/span>, <span style=\"color: #a020f0;\">\"p\"<\/span>)\r\n    imshow(<span style=\"color: #a020f0;\">\"eg.png\"<\/span>)\r\n<span style=\"color: #0000ff;\">end<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_02.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Mask the ship based on the polygon and apply that to the original image.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shipmask = poly2mask(p.Position(:, 1), p.Position(:, 2), height(ship), width(ship));\r\nshipmask = imerode(shipmask, strel(<span style=\"color: #a020f0;\">\"disk\"<\/span>, 2));\r\nfigure\r\ntiledlayout(2, 1, TileSpacing=<span style=\"color: #a020f0;\">\"compact\"<\/span>, Padding=<span style=\"color: #a020f0;\">\"tight\"<\/span>)\r\nnexttile\r\nimshow(shipmask)\r\nnexttile\r\nship = ship.*cast(shipmask, class(ship));\r\nimshow(ship)<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_03.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Get the orientation from the ship mask and apply it to the original image.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shipprops = regionprops(shipmask, <span style=\"color: #a020f0;\">\"Orientation\"<\/span>)\r\nshiphorizontal = imrotate(ship, -shipprops.Orientation);\r\nfigure\r\nimshow(shiphorizontal)<\/pre>\n<pre style=\"font-style: oblique;\">shipprops = \r\n  struct with fields:\r\n\r\n    Orientation: 19.5882\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_04.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Crop to just the ship<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shipbbox = regionprops(logical(shiphorizontal(:, :, 1)), <span style=\"color: #a020f0;\">\"BoundingBox\"<\/span>)\r\nshipthumb = imcrop(shiphorizontal, shipbbox.BoundingBox);\r\nimshow(shipthumb)<\/pre>\n<pre style=\"font-style: oblique;\">shipbbox = \r\n  struct with fields:\r\n\r\n    BoundingBox: [181.5000 320.5000 427 67]\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_05.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<h3>Download Lakeside Satellite Imagery<a name=\"6\"><\/a><\/h3>\n<p>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.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\"><span style=\"color: #0000ff;\">if<\/span> false\r\n<span style=\"color: #228b22;\">% Interactively select limits<\/span>\r\n    w = webmap\r\n    [latlim, lonlim] = wmlimits(w)\r\n<span style=\"color: #0000ff;\">else<\/span>\r\n    latlim = [42.2966508472710 42.3053796253471]\r\n    lonlim = [-71.3800315706306\t-71.3639383165411]\r\n<span style=\"color: #0000ff;\">end<\/span><\/pre>\n<pre style=\"font-style: oblique;\">latlim =\r\n   42.2967   42.3054\r\nlonlim =\r\n  -71.3800  -71.3639\r\n<\/pre>\n<p>Find the ortho layer and download it for this region.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">ortho = wmsfind(<span style=\"color: #a020f0;\">'usgsimageryonly'<\/span>, SearchField=<span style=\"color: #a020f0;\">'serverurl'<\/span>);\r\n[lakeside, lakesidereference] = wmsread(ortho, latlim=latlim, lonlim=lonlim);\r\nfigure\r\ngeoshow(lakeside, lakesidereference)\r\naxis <span style=\"color: #a020f0;\">tight<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_06.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<h3>Prepare the Ship for \"Docking\"<a name=\"8\"><\/a><\/h3>\n<p>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 <img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_eq05820723625559497553.png\" hspace=\"5\" vspace=\"5\" \/> will work well.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shipn75 = imrotate(shipthumb, -75,<span style=\"color: #a020f0;\">'nearest'<\/span>);\r\nfigure\r\nimshow(shipn75);<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_08.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>We'll work in meters and need the dimensions of the ship. A quick web search says the length is <b>399.94m<\/b> and the beam is <b>58.8m<\/b>.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shiplength = 399.94;\r\nshipbeam = 58.8;<\/pre>\n<p>At the same time will check the ratio of length to width with that of our image.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">props = regionprops(logical(shipn75(:,:,1)), <span style=\"color: #a020f0;\">\"MajorAxisLength\"<\/span>, <span style=\"color: #a020f0;\">\"MinorAxisLength\"<\/span>);\r\nimage_shipratio = props.MajorAxisLength\/props.MinorAxisLength\r\nactual_shipratio = shiplength\/shipbeam<\/pre>\n<pre style=\"font-style: oblique;\">image_shipratio =\r\n    6.7436\r\nactual_shipratio =\r\n    6.8017\r\n<\/pre>\n<p>I'm happy with that(!) especially considering I just approximated the ship with a drawn polygon.<\/p>\n<p>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 equations for image dimensions are:<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_eq16977764583318617550.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_eq17791414215949114949.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">imgheightm = shipbeam*sind(15)+shiplength*cosd(15)\r\nimgwidthm = shiplength*sind(15)+shipbeam*cosd(15)<\/pre>\n<pre style=\"font-style: oblique;\">imgheightm =\r\n  401.5309\r\nimgwidthm =\r\n  160.3085\r\n<\/pre>\n<p>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.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">shipoverlay = im2double(shipn75);\r\nshipoverlay(repmat(all(~shipoverlay, 3), [1 1 3])) = NaN;<\/pre>\n<h3>Prepare Lakeside for a Ship<a name=\"13\"><\/a><\/h3>\n<p>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 <b>26986<\/b>.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">p = projcrs(26986)\r\n[latgrid, longrid] = geographicGrid(lakesidereference);\r\n[xgrid, ygrid] = projfwd(p, latgrid, longrid);\r\nfigure\r\nmapshow(xgrid,ygrid,lakeside)\r\naxis <span style=\"color: #a020f0;\">tight<\/span><\/pre>\n<pre style=\"font-style: oblique;\">p = \r\n  projcrs with properties:\r\n\r\n                    Name: \"NAD83 \/ Massachusetts Mainland\"\r\n           GeographicCRS: [1\u00d71 geocrs]\r\n        ProjectionMethod: \"Lambert Conic Conformal (2SP)\"\r\n              LengthUnit: \"meter\"\r\n    ProjectionParameters: [1\u00d71 map.crs.ProjectionParameters]\r\n<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_10.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>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 reproduce this.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">ax = gca;\r\nchart = ax.Children(1);\r\ndatatip(chart, 210495, 894667);<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_09-1.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>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.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">xpt = 210495;\r\nypt = 894667;\r\nmr = maprefcells([xpt xpt+imgwidthm], [ypt-imgheightm ypt], size(shipoverlay, [1 2]), ColumnsStartFrom=<span style=\"color: #a020f0;\">\"north\"<\/span>)<\/pre>\n<pre style=\"font-style: oblique;\">mr = \r\n  MapCellsReference with properties:\r\n\r\n            XWorldLimits: [210495 210655.308527484]\r\n            YWorldLimits: [894265.469065182 894667]\r\n              RasterSize: [432 177]\r\n    RasterInterpretation: 'cells'\r\n        ColumnsStartFrom: 'north'\r\n           RowsStartFrom: 'west'\r\n      CellExtentInWorldX: 0.905697895390393\r\n      CellExtentInWorldY: 0.929469756523354\r\n    RasterExtentInWorldX: 160.3085274841\r\n    RasterExtentInWorldY: 401.530934818089\r\n        XIntrinsicLimits: [0.5 177.5]\r\n        YIntrinsicLimits: [0.5 432.5]\r\n      TransformationType: 'rectilinear'\r\n    CoordinateSystemType: 'planar'\r\n            ProjectedCRS: []\r\n\r\n<\/pre>\n<h3>Stick the Ship!<a name=\"16\"><\/a><\/h3>\n<p>And now we can finally plot the Lakeside image with the Ever Given overlaid on top to scale.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">figure\r\nhold <span style=\"color: #a020f0;\">on<\/span>\r\nmapshow(xgrid, ygrid, lakeside)\r\ns = mapshow(zeros(size(shipoverlay, [1 2])), mr, DisplayType=<span style=\"color: #a020f0;\">\"surface\"<\/span>, CData=shipoverlay);\r\naxis <span style=\"color: #a020f0;\">tight<\/span>\r\nxticks([])\r\nyticks([])<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_10-1.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<\/div>\n","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/pick\/files\/mainEverGivenm_12.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div>\n<p>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.&nbsp;<br \/>\nContents<\/p>\n<p>Get... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/pick\/2021\/04\/02\/get-the-ship-stuck\/\">read more >><\/a><\/p>\n","protected":false},"author":87,"featured_media":12217,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/12187"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/users\/87"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/comments?post=12187"}],"version-history":[{"count":9,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/12187\/revisions"}],"predecessor-version":[{"id":12253,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/12187\/revisions\/12253"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/media\/12217"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/media?parent=12187"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/categories?post=12187"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/tags?post=12187"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}