2D grid Hamiltonian path puzzle solver
1. Review


Code is easier to test (and measure the performance) if it's organized into functions or classes. In this case you have persistent data (W, H, grid, legal, start_pos) that's shared between the various parts of the code, so a class would be the best way to go here.
It would be simpler to write:

EMPTY_SPACE_SYMBOLS = '.'
STARTING_POINT_SYMBOLS = 'Ss'


since Python can iterate over the characters in a string just as easily as over the items in a list.
Assertions should be reserved for detecting programming errors, and not used for reporting problems with the input. That's because assertions can be disabled via the -O option to the Python interpreter, or the PYTHONOPTIMIZE environment variable. For the errors in this code, you should raise exceptions.
It's more general to have a method format_solution that returns a string, rather than a function print_solution that prints it. That way the caller has the option of deciding what to do with it.
Instead of:

[OBSTACLE_SYMBOL for _ in range(self.w)]


write:

[OBSTACLE_SYMBOL] * self.w

Instead of building str_grid from result_grid and then formatting the former, just format result_grid directly:

'\n'.join(' '.join(str(item).ljust(3) for item in row)
          for row in result_grid)

The number "3" only works for small enough puzzles. To make the formatting work in the general case, it should be one more than the maximum number of digits in any number in the path, that is:

len(str(len(path) + 1)) + 1

extend is only called from one place in the code, so it would make sense to inline it at that point.
Coordinates are r and c in some places in the code, and x and y in other places. It would be better to be consistent.


2. Revised code

import time

EMPTY_SPACE_SYMBOLS = '.'
STARTING_POINT_SYMBOLS = 'Ss'
OBSTACLE_SYMBOL = 'X'
DIRS = [(-1, 0), (1, 0), (0, 1), (0, -1)]

class HamiltonSolver:
    """Solver for a Hamilton Path problem."""

    def __init__(self, grid):
        """Initialize the HamiltonSolver instance from a grid, which must be a
        list of strings, one for each row of the grid.

        """
        self.grid = grid
        self.h = h = len(grid)
        self.w = w = len(grid[0])
        if any(len(row) != w for row in grid):
            raise ValueError("Grid is not rectangular")
        self.start = None
        self.legal = set()
        for r, row in enumerate(grid):
            for c, item in enumerate(row):
                if item in STARTING_POINT_SYMBOLS:
                    if self.start is not None:
                        raise ValueError("Multiple starting points")
                    self.start = (r, c)
                elif item in EMPTY_SPACE_SYMBOLS:
                    self.legal.add((r, c))
        if self.start is None:
            raise ValueError("No starting point")

    def format_solution(self, path):
        """Format a path as a string."""
        grid = [[OBSTACLE_SYMBOL] * self.w for _ in range(self.h)]
        for i, (r, c) in enumerate(path, start=1):
            grid[r][c] = i
        w = len(str(len(path) + 1)) + 1
        return '\n'.join(''.join(str(item).ljust(w) for item in row)
                         for row in grid)

    def solve(self):
        """Generate solutions as lists of coordinates."""
        target_path_len = len(self.legal) + 1
        paths = [[self.start]]
        while paths:
            path = paths.pop()
            if len(path) == target_path_len:
                yield path
            r, c = path[-1]
            for dr, dc in DIRS:
                new_coord = r + dr, c + dc
                if new_coord in self.legal and new_coord not in path:
                    paths.append(path + [new_coord])

PUZZLE_GRID = '''
......
.SX...
...X..
.....X
......
XX....
......
'''.split()

def main():
    start_time = time.time()
    n_solutions = 0
    puzzle = HamiltonSolver(PUZZLE_GRID)
    for solution in puzzle.solve():
        if n_solutions == 0:
            print(puzzle.format_solution(solution))
            print("First solution found in {} s"
                  .format(time.time() - start_time))
        n_solutions += 1
    print("{} solution{} found in {} s"
          .format(n_solutions, '' if n_solutions == 1 else 's',
                  time.time() - start_time))

if __name__ == '__main__':
    main()


3. Speeding it up


Each time extend is called, the list of paths under consideration is extended by up to four new paths, each of which is a copy of the current path plus one new grid coordinate. As the search progresses, this results in many paths being created and stored (up to 20 paths in the example from the post).

But in a depth-first search it should be only be necessary to remember a single path at a time, together with the current position in the iteration over the four directions. (This is the "stack of iterators" pattern.)
The code checks to see if the new coordinate is available by calling new_coord in self.legal and new_coord not in path. But path is a list, and testing for membership in a list is not efficient (Python has to compare the new coordinate with every item in the list). For efficient membership testing you need a set.

In this case the code already has a set of legal coordinates. So we can remove a coordinate from this set when we add it to the path, and add it back to the set when we backtrack. This means we only need one membership test, and it also means that the path is complete when the set of legal coordinates is empty.
We can speed up the inner loop by caching method lookups like path.append in local variables.


The revised solve method looks like this:

def solve(self):
    """Generate solutions as lists of coordinates."""
    path = [self.start]
    dirs = [iter(DIRS)]

    # Cache attribute lookups in local variables
    path_append = path.append
    path_pop = path.pop
    legal = self.legal
    legal_add = legal.add
    legal_remove = legal.remove
    dirs_append = dirs.append
    dirs_pop = dirs.pop

    while path:
        r, c = path[-1]
        for dr, dc in dirs[-1]:
            new_coord = r + dr, c + dc
            if new_coord in legal:
                path_append(new_coord)
                legal_remove(new_coord)
                dirs_append(iter(DIRS))
                if not legal:
                    yield path
                break
        else:
            legal_add(path_pop())
            dirs_pop()


On the example from the post this is more than 2½ times as fast as the original code.

Update You asked in comments how this works. Well, the list path contains the current path reached in the search. The list dirs is the same length as path, and each element is an iterator over the list DIRS. Each iterator remembers where we got to in the list of directions at that point in the search, so that when we backtrack we can pop the dirs list and then continue iterating from where we left off.

The control flow is slightly tricky. When the path is extended it reaches the break statement, which breaks out of the for loop (bypassing the else) and so goes back to the top of the while loop, picking up the new coordinates that we just added to the path. But when no further progress can be made, the for loop runs out of directions and enters the else clause. This causes  path and dirs to be popped and then the while loop picks up from the previous point reached in the search.

I hope that's clear. If you're still not quite sure, step through the code in the Python debugger so that you can see exactly what it's doing.