In my game the entire randomly generated world consists of over 1,500,000 tiles. The way that I go about it is a mix of what you said and also what Infinisearch said. The world is split into chunks that hold, in my case, 64x64 tiles each. Only 9 chunks are ever loaded into memory which you can imagine as a 2D array with 3 rows and 3 columns. The chunk in the center of the array is always the chunk that the player is located in. As the player moves throughout the world, old chunks are cached and news chunks are loaded. Chunks contain the information about the individual tiles in that chunk. Tiles are held in a 2D array in the chunk and are represented by integers which correspond to their IDs. When a tile is acted upon, the tile type is found by that ID.

Could you share some pseudo code about how you did that? I'm guessing you use the players world coordinates and either add or subtract from the loaded chunks based on what you need to load

It's a small trick with analytical geometry.

The first part is you need to define what your origin is. So lets say if the origin is set at the center of a chunk at (0,0) local to that chunk. Therfore the origin is also (0,0)

The problem here is that there's not really a way to evenly split off that way, so well have something that looks a bit lopsided. The highest X and Y for each chunk will be 16. The lowest for both will be -15. Notice here that 16 + ( 0 + 1) + |-15| = 32. Note that the zero coordinate counts as a unit, and not nothing.

If you choose an origin point according to that set up, you can then determine how far the origin of each chunk is, then calculate the coordinates for each tile according to the world. I'm sure you can do the rest at this point. I recommend using a sheet of graphing paper to help you determine the equation. I still do that.

And lastly... this is really only for positioning the tiles. When I use a paging system, I typically have everything stored in such a way that it's coordinates are locals of the tile. So when the tile's position is calculated, everything else is calculated based on it's matrix transform and then cached. This may seem slow to some, but honestly you're not calculating these things every frame. And the math is fairly simple.

Lastly... remember that matrix order MATTERS. To make life easier for you.

Scale Matrix first, Rotational Matrix second. Translate Last.

here's how I do it:
a chunk should be at least the size of the player's view
you have 9 chunks loaded at a time (the chunk the player is currently on, plus all 8 neighboring chunks)
your ideal chunk size is such that you can get the 3 new neighboring chunks ready for display in the minimum time it would take the player to cross a single chunk in any direction
if you implement portals or such, this gets much more complicated, unless you think a loading wait is appropriate in the case that the player steps through a portal.

alternative way:
Arrange chunks into "zones" based on how easy it is for the player to move between chunks (or tileset, if you're doing tilesets).
Keep the entire zone ready for display.
Begin loading neighboring zone when player enters a chunk neighboring a chunk with access to neighboring zone.
Chunk and zone sizes are such that an entire zone can be loaded in the time it takes a player to cross a chunk.