How to read .npy file in Mathematica?

I would suggest reading it with numpy and not Mathematica. It seems to me that this is not an exchange format. It is a format meant to be used only by numpy.

Mathematica 12.0 has significantly improved ExternalEvaluate, and now you can transfer data from Python to Mathematica quite efficiently.

ExternalEvaluate["Python",
 "import numpy as np
 x = np.arange(60)
 y=x.reshape(3,4,5)
 y"
]

Normal[%]
(* {{{0, 1, 2, 3, 4}, {5, 6, 7, 8, 9}, {10, 11, 12, 13, 14}, {15,
    16, 17, 18, 19}}, {{20, 21, 22, 23, 24}, {25, 26, 27, 28, 
   29}, {30, 31, 32, 33, 34}, {35, 36, 37, 38, 39}}, {{40, 41, 42, 43,
    44}, {45, 46, 47, 48, 49}, {50, 51, 52, 53, 54}, {55, 56, 57, 58, 
   59}}} *)

This worked on Linux.

str = Import["matrix001.npy", "String", Path -> NotebookDirectory[]];
meta = Flatten@StringCases[str, "{" ~~ __ ~~ "}"];

dims = Flatten@StringCases[meta,
    "(" ~~ z__ ~~ ")" :> ToExpression["{" <> z <> "}"]];
nElems = Times @@ dims;
sizeOfInteger = 8 (* bytes *);

binData = StringTake[str, -sizeOfInteger*nElems ;;];
bstream = StringToStream[binData];
raw = BinaryReadList[bstream, "Integer64"];

On[Assert]
Assert[Length[raw] == nElems, 
     "Did not read the correct number of integers"]
Close[bstream];

data = ArrayReshape[raw, dims]

The idea is to read the entire file as a string. The beginning of the string is text metadata, including the array dimensions. The binary data follows. The code parses out the dimensions, but assumes the data is 8-byte integers, but that could have been parsed out also. Then the code reads the end of the string as a binary stream and reshapes the list into the original shape.

It's completely experimental, but seems to work in this case.