diff --git a/tests/basic_map_operations_test.py b/tests/basic_map_operations_test.py
index 6e216a13..8de981a8 100644
--- a/tests/basic_map_operations_test.py
+++ b/tests/basic_map_operations_test.py
@@ -1,6 +1,6 @@
-import random
+import numpy
 import unittest
-from worldengine.basic_map_operations import distance, index_of_nearest, random_point
+from worldengine.basic_map_operations import distance, index_of_nearest
 
 
 class TestBasicMapOperations(unittest.TestCase):
@@ -8,16 +8,6 @@ class TestBasicMapOperations(unittest.TestCase):
     def setUp(self):
         pass
 
-    def test_random_point(self):
-        for seed in [0, 1, 27, 29, 1939, 1982, 2015]:
-            random.seed(seed)
-            for n in range(10):
-                x, y = random_point(100, 200)
-                self.assertTrue(x >= 0,  "x is within boundaries")
-                self.assertTrue(x < 100, "x is within boundaries")
-                self.assertTrue(y >= 0,  "y is within boundaries")
-                self.assertTrue(y < 200, "y is within boundaries")
-
     def test_distance(self):
         self.assertAlmostEquals(22.360679774997898, distance((0, 0), (10, 20)))
         self.assertAlmostEquals(22.360679774997898, distance((-1, -1), (9, 19)))
diff --git a/worldengine/basic_map_operations.py b/worldengine/basic_map_operations.py
index 7c7b107b..38bfdb5d 100644
--- a/worldengine/basic_map_operations.py
+++ b/worldengine/basic_map_operations.py
@@ -1,9 +1,4 @@
 import math
-import random
-
-
-def random_point(width, height):
-    return random.randrange(0, width), random.randrange(0, height)
 
 
 def distance(pa, pb):
diff --git a/worldengine/cli/main.py b/worldengine/cli/main.py
index 027d56ca..29207c7a 100644
--- a/worldengine/cli/main.py
+++ b/worldengine/cli/main.py
@@ -3,7 +3,6 @@
 import os
 import numpy
 import pickle
-import random
 import worldengine.generation as geo
 from worldengine.common import array_to_matrix, set_verbose, print_verbose
 from worldengine.draw import draw_ancientmap_on_file, draw_biome_on_file, draw_ocean_on_file, \
@@ -390,11 +389,13 @@ def main():
         if not os.path.exists(args.FILE):
             usage("The specified world file does not exist")
 
+    maxseed = 65535  # there is a hard limit somewhere so seeds outside the uint16 range are considered unsafe
     if args.seed is not None:
         seed = int(args.seed)
+        assert 0 <= seed <= maxseed, "Seed has to be in the range between 0 and %s, borders included." % maxseed
     else:
-        seed = random.randint(0, 65535)#RNG initialization is done automatically
-    random.seed(seed)
+        seed = numpy.random.randint(0, maxseed)  # first-time RNG initialization is done automatically
+    numpy.random.seed(seed)
 
     if args.world_name:
         world_name = args.world_name
diff --git a/worldengine/drawing_functions.py b/worldengine/drawing_functions.py
index bd9f6202..46a1bd71 100644
--- a/worldengine/drawing_functions.py
+++ b/worldengine/drawing_functions.py
@@ -6,7 +6,6 @@
 
 import math
 import numpy
-import random
 import sys
 import time
 from worldengine.common import get_verbose
@@ -133,8 +132,7 @@ def is_inner_border(pos):
 
 
 def _find_mountains_mask(world, factor):
-    _mask = [[False for x in range(factor * world.width)] for y in
-             range(factor * world.height)]
+    _mask = numpy.full((factor * world.height, factor * world.width), False, dtype=object)
     for y in range(factor * world.height):
         for x in range(factor * world.width):
             if world.is_mountain((int(x / factor), int(y / factor))):
@@ -142,7 +140,7 @@ def _find_mountains_mask(world, factor):
                                            radius=3,
                                            predicate=world.is_mountain))
                 if v > 32:
-                    _mask[y][x] = v / 4
+                    _mask[y, x] = v / 4.0  # force conversion to float, Python 2 will *not* do it automatically
     return _mask
 
 
@@ -515,7 +513,7 @@ def _draw_savanna(pixels, x, y):
 
 
 # TODO: complete and enable this one
-def _dynamic_draw_a_mountain(pixels, x, y, w=3, h=3):
+def _dynamic_draw_a_mountain(pixels, rng, x, y, w=3, h=3):
     # mcl = (0, 0, 0, 255)  # TODO: No longer used?
     # mcll = (128, 128, 128, 255)
     mcr = (75, 75, 75, 255)
@@ -530,7 +528,7 @@ def _dynamic_draw_a_mountain(pixels, x, y, w=3, h=3):
         max_leftborder = int(bottomness * w * 1.33)
         if not last_leftborder == None:
             max_leftborder = min(max_leftborder, last_leftborder + 1)
-        leftborder = int(bottomness * w) + random.randint(-2, 2)/2
+        leftborder = int(bottomness * w) + rng.randint(-2, 2)/2
         if leftborder < min_leftborder:
             leftborder = min_leftborder
         if leftborder > max_leftborder:
@@ -557,7 +555,7 @@ def _dynamic_draw_a_mountain(pixels, x, y, w=3, h=3):
         max_modx = int(bottomness * w * 1.33)
         if not last_modx == None:
             max_modx = min(max_modx, last_modx + 1)
-        modx = int(bottomness * w) + random.randint(-2, 2)/2
+        modx = int(bottomness * w) + numpy.random.randint(-2, 2)/2
         if modx < min_modx:
             modx = min_modx
         if modx > max_modx:
@@ -603,7 +601,7 @@ def draw_ancientmap(world, target, resize_factor=1,
                     sea_color=(212, 198, 169, 255),
                     draw_biome = True, draw_rivers = True, draw_mountains = True,
                     draw_outer_land_border = False, verbose=get_verbose()):
-    random.seed(world.seed * 11)  
+    rng = numpy.random.RandomState(world.seed)  # create our own random generator
 
     if verbose:
         start_time = time.time()
@@ -893,7 +891,7 @@ def _anti_alias_point(x, y):
                     if len(world.tiles_around_factor(resize_factor, (x, y),
                                                      radius=r,
                                                      predicate=on_border)) <= 2:
-                        if random.random() <= .5:
+                        if rng.random_sample() <= .5:
                             _draw_temperate_forest1(target, x, y, w=w, h=h)
                         else:
                             _draw_temperate_forest2(target, x, y, w=w, h=h)
diff --git a/worldengine/generation.py b/worldengine/generation.py
index 2453ef17..099f5134 100644
--- a/worldengine/generation.py
+++ b/worldengine/generation.py
@@ -194,31 +194,45 @@ def _around(x, y, width, height):
 def generate_world(w, step):
     if isinstance(step, str):
         step = Step.get_by_name(step)
-    seed = w.seed
 
     if not step.include_precipitations:
         return w
 
+    # Prepare sufficient seeds for the different steps of the generation
+    rng = numpy.random.RandomState(w.seed)  # create a fresh RNG in case the global RNG is compromised (i.e. has been queried an indefinite amount of times before generate_world() was called)
+    sub_seeds = rng.randint(0, 4294967295, size=100)  # sys.maxsize didn't quite work
+    seed_dict = {
+                 'PrecipitationSimulation': sub_seeds[ 0],  # after 0.19.0 do not ever switch out the seeds here to maximize seed-compatibility
+                 'ErosionSimulation':       sub_seeds[ 1],
+                 'WatermapSimulation':      sub_seeds[ 2],
+                 'IrrigationSimulation':    sub_seeds[ 3],
+                 'TemperatureSimulation':   sub_seeds[ 4],
+                 'HumiditySimulation':      sub_seeds[ 5],
+                 'PermeabilitySimulation':  sub_seeds[ 6],
+                 'BiomeSimulation':         sub_seeds[ 7],
+                 '':                        sub_seeds[99]
+    }
+
     # Precipitation with thresholds
-    PrecipitationSimulation().execute(w, seed)
+    PrecipitationSimulation().execute(w, seed_dict['PrecipitationSimulation'])
 
     if not step.include_erosion:
         return w
-    ErosionSimulation().execute(w, seed)
+    ErosionSimulation().execute(w, seed_dict['ErosionSimulation'])  # seed not currently used
     if get_verbose():
         print("...erosion calculated")
 
-    WatermapSimulation().execute(w, seed)
+    WatermapSimulation().execute(w, seed_dict['WatermapSimulation'])  # seed not currently used
 
     # FIXME: create setters
-    IrrigationSimulation().execute(w, seed)
-    TemperatureSimulation().execute(w, seed)
-    HumiditySimulation().execute(w, seed)
+    IrrigationSimulation().execute(w, seed_dict['IrrigationSimulation'])  # seed not currently used
+    TemperatureSimulation().execute(w, seed_dict['TemperatureSimulation'])
+    HumiditySimulation().execute(w, seed_dict['HumiditySimulation'])  # seed not currently used
 
     
-    PermeabilitySimulation().execute(w, seed)
+    PermeabilitySimulation().execute(w, seed_dict['PermeabilitySimulation'])
 
-    cm, biome_cm = BiomeSimulation().execute(w, seed)
+    cm, biome_cm = BiomeSimulation().execute(w, seed_dict['BiomeSimulation'])  # seed not currently used
     for cl in cm.keys():
         count = cm[cl]
         if get_verbose():
diff --git a/worldengine/plates.py b/worldengine/plates.py
index 05430a89..9978ccbd 100644
--- a/worldengine/plates.py
+++ b/worldengine/plates.py
@@ -2,7 +2,6 @@
 # extension which is not available when using this project from jython
 
 import platec
-import random
 import time
 import numpy
 
@@ -64,7 +63,7 @@ def world_gen(name, width, height, seed, num_plates=10, ocean_level=1.0,
 
     if verbose:
         start_time = time.time()
-    add_noise_to_elevation(world, random.randint(0, 4096))
+    add_noise_to_elevation(world, numpy.random.randint(0, 4096))  # uses the global RNG; this is the very first call to said RNG - should that change, this needs to be taken care of
     if verbose:
         elapsed_time = time.time() - start_time
         print("...plates.world_gen: elevation noise added. Elapsed time " +
diff --git a/worldengine/simulations/erosion.py b/worldengine/simulations/erosion.py
index f73d7a20..dbac06a6 100644
--- a/worldengine/simulations/erosion.py
+++ b/worldengine/simulations/erosion.py
@@ -32,6 +32,7 @@ def _numpy_to_matrix(numpy_array):
 
     This is used because currently we do not know how to serialize numpy
     arrays :( with pickle. In the future we will use pytables/hdf5"""
+    # TODO: Is this still true? Pickle didn't seem to cause problems for me. /tcld
 
     width = numpy_array.shape[0]
     height = numpy_array.shape[1]
diff --git a/worldengine/simulations/hydrology.py b/worldengine/simulations/hydrology.py
index 3f1ecc74..e0424075 100644
--- a/worldengine/simulations/hydrology.py
+++ b/worldengine/simulations/hydrology.py
@@ -57,8 +57,8 @@ def droplet(world, pos, q, _watermap):
 
         _watermap_data = numpy.zeros((world.height, world.width), dtype=float)
         for i in range(n):
-            x, y = world.random_land()
-            if True and world.precipitation['data'][y, x] > 0:
+            x, y = world.random_land()  # will return None for x and y if no land exists
+            if x is not None and world.precipitation['data'][y, x] > 0:
                 droplet(world, (x, y), world.precipitation['data'][y, x],
                         _watermap_data)
         _watermap = dict()
diff --git a/worldengine/simulations/irrigation.py b/worldengine/simulations/irrigation.py
index 6c8bc4a1..5a588920 100644
--- a/worldengine/simulations/irrigation.py
+++ b/worldengine/simulations/irrigation.py
@@ -5,7 +5,7 @@ class IrrigationSimulation(object):
     def is_applicable(world):
         return world.has_watermap() and (not world.has_irrigation())
 
-    def execute(self, world, seed):#seed is currently not used
+    def execute(self, world, seed):
         world.irrigation = self._calculate(world)
 
     @staticmethod
diff --git a/worldengine/simulations/permeability.py b/worldengine/simulations/permeability.py
index 8fbae206..6ed21756 100644
--- a/worldengine/simulations/permeability.py
+++ b/worldengine/simulations/permeability.py
@@ -1,6 +1,5 @@
 from worldengine.simulations.basic import find_threshold_f
 from noise import snoise2
-import random
 import numpy
 
 
@@ -21,8 +20,9 @@ def execute(self, world, seed):
 
     @staticmethod
     def _calculate(seed, width, height):
-        random.seed(seed * 37)
-        base = random.randint(0, 4096)
+        rng = numpy.random.RandomState(seed)  # create our own random generator
+        base = rng.randint(0, 4096)
+
         perm = numpy.zeros((height, width), dtype=float)
 
         octaves = 6
diff --git a/worldengine/simulations/precipitation.py b/worldengine/simulations/precipitation.py
index 77215b01..a6d48136 100644
--- a/worldengine/simulations/precipitation.py
+++ b/worldengine/simulations/precipitation.py
@@ -1,4 +1,3 @@
-import random
 import time
 import numpy
 from noise import snoise2
@@ -32,9 +31,10 @@ def execute(self, world, seed):
     @staticmethod
     def _calculate(seed, width, height):
         """Precipitation is a value in [-1,1]"""
+        rng = numpy.random.RandomState(seed)  # create our own random generator
+        base = rng.randint(0, 4096)
+
         border = width / 4
-        random.seed(seed * 13)
-        base = random.randint(0, 4096)
         precipitations = numpy.zeros((height, width), dtype=float)
 
         octaves = 6
diff --git a/worldengine/simulations/temperature.py b/worldengine/simulations/temperature.py
index e2594787..84d998aa 100644
--- a/worldengine/simulations/temperature.py
+++ b/worldengine/simulations/temperature.py
@@ -1,5 +1,4 @@
 from worldengine.simulations.basic import find_threshold_f
-import random
 import numpy
 
 
@@ -31,8 +30,8 @@ def _calculate(world, seed, elevation, mountain_level):
         width = world.width
         height = world.height
 
-        random.seed(seed * 7)
-        base = random.randint(0, 4096)
+        rng = numpy.random.RandomState(seed)  # create our own random generator
+        base = rng.randint(0, 4096)
         temp = numpy.zeros((height, width), dtype=float)
 
         from noise import snoise2
diff --git a/worldengine/world.py b/worldengine/world.py
index 7110a925..9a687f6e 100644
--- a/worldengine/world.py
+++ b/worldengine/world.py
@@ -12,7 +12,6 @@
     WarmTemperateWetForest, TropicalDesert, TropicalDesertScrub, TropicalDryForest, \
     TropicalMoistForest, TropicalRainForest, TropicalThornWoodland, TropicalWetForest, \
     TropicalVeryDryForest, biome_index_to_name, biome_name_to_index
-from worldengine.basic_map_operations import random_point
 import worldengine.protobuf.World_pb2 as Protobuf
 from worldengine.step import Step
 from worldengine.common import _equal
@@ -332,11 +331,12 @@ def contains(self, pos):
     #
 
     def random_land(self):
-        x, y = random_point(self.width, self.height)
-        if self.ocean[y, x]:#TODO: this method should get a safer/quicker way of finding land!
-            return self.random_land()
-        else:
-            return x, y
+        if self.ocean.all():
+            return None, None  # return invalid indices if there is no land at all
+        lands = numpy.invert(self.ocean)
+        lands = numpy.transpose(lands.nonzero())  # returns a list of tuples/indices with land positions
+        y, x = lands[numpy.random.randint(0, len(lands))]  # uses global RNG
+        return x, y
 
     def is_land(self, pos):
         return not self.ocean[pos[1], pos[0]]#faster than reversing pos or transposing ocean