I have a pretty involved solution that requires you specify some hyperparameters of your model. Not being an expert in the relation between snake length and age, so I picked some values that seemed vaguely possible. But, by all means feel to tweak them. This series of functions is adapted and expanded from Gelman and Hill (2007).

First, we need a function to simulate data sets with the properties you want. Here is where I "invented" some numbers (under the comment #Specifying hyperparameters. Feel free to adjust accordingly based on your knowledge base and expectations.

#Creating fake dataframe for simulations 
fake<-function(J, K){
  snake<-(1:(J*K))          #creates a snake or case id 
  habitat<-rep(1:J, K)      #creates a balanced grouping variable   
  age<-rnorm(J*K, 2.5, .75) #creates a random normal age variable 
  #Specifying hyperparameters 
  mu.a.true<-12 #grand intercept 
  mu.g1.true<-3 #grand age slope 
  sigma.y.true<-1.65        #residual sd 
  sigma.a.true<-2.30        #intercept sd
  sigma.g1.true<-2.15       #age slope sd 
  #habitat-level parameters 
  a.true<-rnorm(J, mu.a.true, sigma.a.true) #vector of random intercepts
  g1.true<-rnorm(J, mu.g1.true, sigma.g1.true) #vector of random slopes  
  #dependent variable 
  y<-rnorm(J*K, a.true[habitat]+g1.true[habitat]*age, sigma.y.true)
  return(data.frame(snake, y, habitat, age))
}

and if you run print(fake(5,5)) you get the following simulated data (where y = length):

> print(fake(5,5))
   snake        y habitat       age
1      1 21.85790       1 1.9030948
2      2 27.73040       2 2.4382916
3      3 11.60743       3 0.7640231
4      4 15.63832       4 2.9965471
5      5 26.83866       5 1.8319489
6      6 24.65482       1 3.0666732
7      7 27.42804       2 2.2352739
8      8 16.98225       3 3.5354592
9      9 15.32236       4 3.0284722
10    10 25.22085       5 2.4662031
11    11 18.39988       1 2.2474545
12    12 21.85225       2 1.2652154
13    13 14.12972       3 2.6137810
14    14 17.69987       4 1.3900062
15    15 33.13150       5 3.2806024
16    16 28.80915       1 4.6216763
17    17 35.61110       2 2.5738612
18    18 15.78296       3 2.5815051
19    19 17.76041       4 3.2101400
20    20 29.23588       5 2.5638092
21    21 21.51020       1 2.7942239
22    22 32.10491       2 2.6981291
23    23 10.05619       3 1.7800981
24    24 16.96727       4 3.0558285
25    25 34.61646       5 3.3928710

Up next? We need to simulate a bunch of data and see how many data sets "recover" a significant effect for age given the hyperparameters we input in our data simulation function. Again, those hyperparameters can be adjusted to better fit your expectations.

mixed.power<-function(J, K, n.sims=1000){
  signif<-rep(NA, n.sims) #note that you can specify number of simulations - default is 1000
  pb<-winProgressBar(title="Progress", min=0, max=100, width=300)    #if you want to watch the progress
  require(lme4)
  require(arm)
  for(s in 1:n.sims){
    fake.data<-fake(J, K)                                               #calls in data simulation function
    lme.power<-lmer(y~1+age+(1+age|habitat), data=fake.data)            #estimates mixed effect model using each simulated dataset
    theta.hat<-fixef(lme.power)["age"]                                  #saves age coefficients from each simulated dataset 
    theta.se<-se.fixef(lme.power)["age"]                                #saves standard error of age coefficients from each simulated dataset
    signif[s]<-ifelse((abs(theta.hat)-qt(.975, J*K-1)*theta.se)>0, 1, 0)#assigns value of 1 to significant coefficients 0 to ns coefficients
    setWinProgressBar(pb, s/n.sims*100, title=paste(round(s/n.sims*100, 0), "% done"))
  }
  close(pb)
  power<-mean(signif, na.rm=T)#calculates proportion of significant models out of # of simulated datasets... 
  return(power)
}

You may get some warnings with these codes when models don't converge in your simulations. I find that happens more often when dealing with small samples like the ones you are proposing. (I tend to ignore these convergence warnings)

You can pause here and calculate power for a given number of habitats (J) with a given number of snakes per habitat (K). In that case you can run mixed.power(5,10) and see the following returned:

> mixed.power(5, 10)
[1] 0.804

Or perhaps more instructively you can go all out and create a power plot to get a sense of how many snakes per habitat you should consider gathering. Here is the function for that:

graph.power<-function(J, max.K){
  require(ggplot2)
  KK<-seq(3, max.K, by=1)       #will ieterate from 3 cases to the max.K specified in the function above
  Y<-rep(NA, length(KK))        #Empty vector to contain power estimates from mixed.power () function
  JJ<-rep(J, length(KK))        #Repeates the number of habitats you specified to equal the length of the KK vector 
  pb2<-winProgressBar(title="Overall Progress", min=0, max=100, width=300)  #for an overall progress bar
  for(i in 1:length(KK)){
    Y[i]<-mixed.power(JJ[i], KK[i]) #runs mixed power function and stores results one at a time
    setWinProgressBar(pb2, i/length(KK)*100, title=paste(round(i/length(KK)*100, 0), "% Overall complete"))
  }
  close(pb2)
  DF<-as.data.frame(cbind(KK, Y)) #creates data frame for plotting
  colnames(DF)<-c("KK", "Y")
  #Code below provides a basic power plot with a horizontal line at .80 
  g1<-ggplot(aes(x=KK, y=Y), data=DF) 
  g2<-g1+geom_smooth(se=F)+geom_hline(yintercept = .8, lty="dashed", col="red", lwd=2)
  g3<-g2+xlab("Sample Size per Habitat")+ylab("Power")
  g4<-g3+ggtitle("Power Analysis")
  return(g4)
}

When you use this function, the only difference is that now the second number is the maximum number of snakes you want to consider examining in a power analysis. So in this example, let's say I wanted to look at a power curve that considered up to 15 snakes per habitat (holding the number of habitats constant at 5). Here is what that would look like: graph.power(5,15)... which would return the following image:

Now you have the ability to assess the sample size at which the power curve crosses .80 (~9 snakes per 5 habitats in this case).

This is NOT the only way to approach a power analysis in a mixed model, but I tend to like it as it gives me a great deal of control over the various hyperparameters and allows me to plot various power curves to get a sense of how changing those hyperparameters in addition to sample and cluster size change the power of my model.

A note: if you want to play around to get a sense of how the functions work or to try to identify some good starting values you may want to lower the number of simulations in the mixed.power() function to reduce computation time prior to running your "official" power analyses with the larger set of simulated data sets.