You can filter using contour area then apply morpholgical closing to fill the small holes in the image. Here's the result:

import cv2

# Load image, convert to grayscale, Gaussian blur, Otsu's threshold
image = cv2.imread('1.jpg')
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
blur = cv2.GaussianBlur(gray, (3,3), 0)
thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)[1]

# Filter using contour area and remove small noise
cnts = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
cnts = cnts[0] if len(cnts) == 2 else cnts[1]
for c in cnts:
    area = cv2.contourArea(c)
    if area < 5500:
        cv2.drawContours(thresh, [c], -1, (0,0,0), -1)

# Morph close and invert image
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5,5))
close = 255 - cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel, iterations=2)

cv2.imshow('thresh', thresh)
cv2.imshow('close', close)
cv2.waitKey()

import numpy as np
import matplotlib.pyplot as plt

# image posted by OP
URL = "https://i.sstatic.net/Pa7Io.png"

# Read image
from skimage import io
from skimage.filters import threshold_otsu
from skimage.color import rgb2gray

image = rgb2gray(io.imread(URL)[21:899, 555:1125, :3]) #index cut the needed part

# Convert to True-False
_2color = image > threshold_otsu(image)

# accumulated true from left to right. Is true if any pixel to the left is True
trueLeft = np.logical_or.accumulate(_2color, axis=1)
# accumulated true from right to left. Is True if any pixel to the right is True
trueRight = np.logical_or.accumulate(_2color[:, ::-1], axis=1)[:, ::-1]

# True if has any true in the left and right
orImage = trueLeft * trueRight

plt.imshow(orImage)
plt.show()

Use cv2.connectedComponentsWithStats (doc):

import cv2
import numpy as np

# Start by finding all of the connected components (white blobs in your image).
# 'im' needs to be grayscale and 8bit.
nb_blobs, im_with_separated_blobs, stats, _ = cv2.connectedComponentsWithStats(im)
# im_with_separated_blobs is an image where each detected blob has a different pixel value ranging from 1 to nb_blobs - 1.
# The background pixels have value 0.

im_with_separated_blobs looks like this :

# 'stats' (and the silenced output 'centroids') provides information about the blobs. See the docs for more information. 
# Here, we're interested only in the size of the blobs :
sizes = stats[:, cv2.CC_STAT_AREA]
# You can also directly index the column with '-1' instead of 'cv2.CC_STAT_AREA' as it's the last column.

# A small gotcha is that the background is considered as a blob, and so its stats are included in the stats vector at position 0.

# minimum size of particles we want to keep (number of pixels).
# here, it's a fixed value, but you can set it as you want, eg the mean of the sizes or whatever.
min_size = 150  

# create empty output image with will contain only the biggest composents
im_result = np.zeros_like(im_with_separated_blobs)

# for every component in the image, keep it only if it's above min_size.
# we start at 1 to avoid considering the background
for index_blob in range(1, nb_blobs):
    if sizes[index_blob] >= min_size:
        im_result[im_with_separated_blobs == index_blob] = 255

im_result looks like this :

Bonus : An alternative way of coding the same process, taking advantage of indexing a numpy array with another array and of the numpy.where function. It's probably more Pythonic insofar that it avoids the for loop, but is also possibly less flexible rearding other filtering conditions (and perhaps less readable) :

_, im_with_separated_blobs, stats, _ = cv2.connectedComponentsWithStats(im)

sizes = stats[:, cv2.CC_STAT_AREA]

im_result = np.where(sizes[im_with_separated_blobs] >= min_size, im, 0)

The key here is to rely on the indexing of a numpy array with a 2D array. sizes[im_with_separated_blobs] is a 2D numpy array with the same shape as im_with_separated_blobs, but where every value has been replaced by sizes[value], i.e. the size of that blob.

This might be a good contribution to the scikit-image library itself, but for now, you need to roll your own. As suggested by Christoph, you can subtract the result of remove_small_objects from the original image to remove large objects. So, something like:

def filter_objects_by_size(label_image, min_size=0, max_size=None):
    small_removed = remove_small_objects(label_image, min_size)
    if max_size is not None:
        mid_removed = remove_small_objects(small_removed, max_size)
        large_removed = small_removed - mid_removed
        return large_removed
    else:
        return small_removed

You can try using erosion/dilation from OpenCV. Here is a quick example, edit the parameters as needed.

import cv2
import numpy as np

img = cv2.imread('img.png')
blurred_img = cv2.medianBlur(img, 5)

kernel = np.ones((3,3),np.uint8)
erosion = cv2.erode(blurred_img, kernel, iterations=1)
output = cv2.dilate(erosion, kernel, iterations=1)

cv2.imwrite('output.png', output)