Although there are >3,700 species of snakes worldwide (Uetz et al., 2020), 600–800 of which are medically important (Uetz et al., 2020), we chose to initially focus on 45 of the most well-represented species, each with ≥500 photos per species (i.e., per class). We gathered a total of 82,601 images from open online citizen science biodiversity platforms (iNaturalist, HerpMapper) and photo sharing sites (Flickr). The photos were labeled by users of these platforms; either by the user who uploaded the photo (Flickr, HerpMapper) or through a consensus reached by users who viewed the photo and submitted an identification (iNaturalist; see Hochmair et al., 2020 for a synopsis). The data from iNaturalist were collected Oct 24, 2018 using the iNaturalist export tool (https://www.inaturalist.org/observations/export) with the parameters quality_grade=research&identifications=any&captive=false&taxon_id=85553. The HerpMapper data were provided via a partner request on 19 Sept 2018 (see https:// www.herpmapper.org/data-request). The Flickr data were collected Nov 5, 2018 using a python script (https://github.com/cam4ani/snakes/blob/master/get_flickr_data.ipynb). The full training dataset can be accessed at (https://datasets.aicrowd.com/aws-eu-central1/aicrowd-static/datasets/snake-species-identification-challenge/train.tar.gz)

In order to improve the accuracy of labels used for image classification, we employed several best practices. For the iNaturalist data, only “research grade” observations, which require a community-supported and agreed-upon identification at the species taxonomic rank, were used. Additionally, we selected a subset (N = 336) of these images for additional label validation by human experts (Durso et al., 2021) and found that just five (1.5%) were misidentified (these have been corrected on iNaturalist). HerpMapper is used primarily by experienced enthusiasts with a lot of experience in snake identification; we found no misidentifications in the subset (N = 200) we examined. Finally, we used only scientific names in our Flickr search to target photographers with a more serious interest in biodiversity. We found no misidentifications in the subset (N = 63) we examined, although one image showed only the habitat without an actual snake. Subset sample sizes were chosen to represent 0.5% of the total dataset from each source. Human experts (N = 250) were recruited via social media, targeting Facebook groups that specialize in snake identification (Durso et al., 2021) as well as by email or private messages over iNaturalist to top identifiers. Because of privacy issues, we were unable to collect demographic information about this community of experts. Their accuracy in aggregate was 43% at the species level, 56% at the genus level, and 75% at the family level, with important variation among individuals, snake species, and global region (Durso et al., 2021).

The distribution of photos among classes was unequal: the class with the most photos had 11,092, while the class with the fewest photos had 517. Some species classes have high intraclass variance (Figure 1), due to geographic and ontogenetic variation (e.g., Manier, 2004), and to color and pattern polymorphism (e.g., Shannon and Humphery, 1963; Cox et al., 2012), and others have very low interclass variance (Figure 1), due to morphological similarity among closely related species as well as inter-species and even inter-family mimicry (e.g., Sweet, 1985; Akcali and Pfennig, 2017). The 45 species we used are found largely in North America (42 species in 19 genera; see Appendix I for a full list), with two Eurasian species (Hierophis viridiflavus and Natrix natrix) and one Central/South American species (Boa constrictor). An interactive version of Figure 1 is available at https://chart-studio.plotly.com/∼amdurso/1/#/.

We used the platform AICrowd, which takes a collaborative approach to the development of algorithms, by inviting data science experts and enthusiasts to collaboratively solve real-world problems by participating in challenges in which the solutions are automatically evaluated in real-time. On January 21, 2019, we launched a “Snake Species Identification Challenge”. The first round lasted until May 31, 2019, and the second round (during which live code collection was implemented) until July 31, 2019. We offered prizes as incentives for the best algorithms submitted (a travel grant and co-authorship on this manuscript). A total of 24 participants made 356 submissions, resulting in five algorithms with an F1 ≥ 0.75 and a top score of F1 = 0.861 with a log-loss = 0.53.

We evaluated the identification accuracy of the submitted algorithms using two test datasets (TD2 and TD3), both distinct from the training dataset described above and its subset (TD1) used for pre-submission testing by the challenge participant (Table 1). One (TD2) was made up of 42,688 photos of the same 45 species of snakes submitted to iNaturalist between January 1, 2019 and September 2, 2019 (after the beginning of our challenge; range 23–6,228 photos per class). The other (TD3) was made up of 248 undisclosed images from 27 classes (1–10 images per class) that were collected from private individuals. Many of the images in this dataset were purposefully chosen to be as difficult to identify as possible—e.g., low resolution, out of focus, with the snake filling only a small part of the frame, and/or obscured by vegetation. The identity of species in these images was confirmed by a herpetologist (A. Durso) and they were used in a citizen science challenge where they were presented to participants recruited from online snake identification communities (largely Facebook snake identification groups and iNaturalist) who suggested species identifications, resulting in 68–157 labels of 5–47 classes per image (Durso et al., 2021). We further subdivided TD3 to take all 45 classes into account (TD3a; for fairer comparison with human labeling, because humans were allowed to choose any of the >3,700 snake species classes from a list), and to evaluate only the 27 classes in common (TD3b; for fairer comparison with TD2).

Applying large, deep convolutional neural networks for image classification is a well-studied problem (Krizhevsky et al., 2012; Huang et al., 2017; He et al., 2019). The top algorithm made use of incremental learning in neural networks and incorporated elements of EfficientNet from Google Brain (Tan and Le, 2019), a pre-trained network from ILSRVC (Russakovsky et al., 2015), discriminative learning, cyclic learning rates and automated image object detection (Figure 2).

As shown by Kornblith et al. (2019), models trained on a standard data distribution generalize better than the models trained from scratch. The ImageNet Large Scale Visual Recognition Competition (ILSVRC) dataset is the most widely used dataset for benchmarking image classifiers, comprising 1.2 million images classified into 1,000 different classes. The winning solution applied incremental learning on a pre-trained EfficientNet network. Specifically, this involved retaining what the model has learned from the ILSVRC dataset and performing incremental learning on the snake species domain. The final layer from the pretrained network was removed and replaced with a domain specific head, a fully connected layer of size 45, each representing the probability of the snakes being in a particular class. Discriminative learning strategy was also used to train the network. Specifically, different layers of the network are responsible for capturing different types of information (Yosinski et al., 2014) and discriminative learning allows us to set the rate at which different components of the network learn. The initial layers are trained at much lower learning rates to inhibit the loss of learned information while the final layers are trained at higher learning rates. A general update to model parameters ⍵ at time step t looks like: ω t =   ω t − 1 − λ . Δ ω J ( ω ) where λ denotes the learning rate and ΔωJ(ω) denotes the gradient with respect to the model’s objective function. In discriminative learning, the model is split into N components, {ω ¹, ω ²,…ω ^N} where ω ⁿ contains the layers of the nth component of the model. Each component can have any number of layers. An update to model parameters then becomes: ω t n =   ω t − 1 n − λ n . Δ ω J ( ω ) where λ ⁿ denotes the learning rate of the nth component of the model.

In the first attempt, an image classifier was trained using progressive resizing, starting with Densenet121 (Huang et al., 2017) architecture with a modified focal loss function (Lin et al., 2017) for multi-class image classification and resizing image sizes from (256,256) → (384,384) → (512,512) → (768,768) → (1024,1024). Using this as an initial method, the scores plateaued at F1 ∼0.67. As an alternative, a new image classifier was trained from scratch using Resnet152 with modifications as suggested by Bag of tricks (He et al., 2019) called XResnet152 (152 indicating the number of layers). This time, the scores plateaued at ∼0.75.

Adding a preprocessing pipeline to predict the four coordinates of the corners of the box bounding the snake itself (which may be any size and shape and in any orientation within the image, against any background) and crop the images, as well as handling orientation variance, was accomplished by annotating 30–32 images from the training dataset from each species category using an annotation tool called sloth (https://github.com/cvhciKIT/sloth/). Pipelining preprocessing and training with XResnet architecture together increases the accuracy to 0.78–0.79.

Finally, a pre-trained EfficientNet (Tan and Le, 2019), which balances the depth of the architecture, the width of the architecture [Mobile Inverted Bottleneck Convolutional block (Sandler et al., 2018) with Swish activation function (Ramachandran et al., 2017)], image resolution and uses appropriate drop-outs, was fine-tuned using preprocessed images. Our winning algorithm carefully tuned hyper-parameters for EfficientNetB0 (Supplementary Figure S1) and used the exact same parameters for EfficientNetB5 (Supplementary Figure S2), although it is expected that higher accuracy is possible using the more time-intensive checkpoints/pretrained weights for B6 and B7 and by fine-tuning EfficientNetB6/B7 pre-trained on ImageNet. Some of the best hyperparameters required to train the EfficientNet include: 1) taking off the final layer of EfficientNetB5 and adding a single layer of size 45 (fastai adds an additional layer, which was not optimal for this case); 2) using the LabelSmoothingCrossEntropy Loss function; 3) using RMSProp optimizer with centered = True; 4) setting bn_wd = false, no batchnorm weight decay; 5) training with discriminative learning; 6) using resize method as SQUISH and not CROP; 7) no mixup augmentation; 8) use rotation augmentation (rotating the training image from −90 to +90 with a probability of 0.8) as shown in Figure 3; and 9) training on 95% of the dataset instead of the 80–20 split.

The network was trained on a single Tesla V100 GPU with a batch size of 8 for 12 epochs. The learning rate schedule is shown in Figure 4. Each epoch took approximately 55 min for completion. Other approaches that were tried but did not work include building a genus classifier to predict the snake genus and then the snake species within genus, focal loss, a weighted sampler/oversampling and batch accumulation (updates happen every n epochs, where n > 1).

The algorithm and readme are available at https://github.com/GokulEpiphany/contests-final-code/tree/master/aicrowd-snake-species.

A confusion matrix generated using TD1 (a withheld subset of the training data) is shown in Figure 5. Species classification accuracy ranged from 96% for Rhinocheilus lecontei to 35% for Pantherophis spiloides. Interclass confusion was 0 for 1332/2025 (66%) of class pairs. F1 for this dataset was 0.83 and log-loss was 0.49 (Table 1). Relatively high confusion remains between some similar species pairs. In particular, the three putative species of the Pantherophis obsoletus (North American ratsnake) complex (Burbrink, 2001) were frequently confused (Table 2).

Among species that are uncontroversially delimited, there were seven species pairs with confusion >10% in TD1 (max = 26% for Thamnophis elegans as Thamnophis sirtalis; Table 3). The highest confusion between species in different genera was Natrix natrix as Thamnophis sirtalis (7%). Although these genera are found in different hemispheres, confusion among co-occurring genera exceeded 5% for Coluber constrictor as Pantherophis alleghaniensis (6%) and for Lichanura trivirgata as Crotalus atrox (5%). The last is particularly troubling because Lichanura trivirgata is a harmless boid whereas Crotalus atrox is a large, potentially dangerous rattlesnake; these species co-occur in the Sonoran Desert and are probably frequently photographed against similar backgrounds (Rorabaugh, 2008). Other frequently confused intergeneric and intercontinental species pairs were Hierophis viridiflavus (as Thamnophis sirtalis and as Masticophis flagellum; both 5%) and Natrix natrix (as Nerodia rhombifer; 5%).

In our second test dataset (TD2, containing images for all 45 classes taken from iNaturalist after the challenge had started), F1 = 0.83 and log-loss = 0.66 (Table 1). Species classification accuracy ranged from 97% for Pantherophis guttatus to 61% for Pantherophis alleghaniensis. Interclass confusion was 0 for 1118/2025 (55%) of class pairs. Again, the most frequently confused were members of the Pantherophis obsoletus complex (Table 2).

Among species that are uncontroversially delimited, there were five species pairs with confusion >10% in TD2 (Table 3). The highest confusion among species in different genera was Pituophis catenifer as Lichanura trivirgata (14%). As in TD1, Pituophis catenifer and Lichanura trivirgata co-occur in the Sonoran Desert and are probably frequently photographed against similar backgrounds, although they are dissimilar in appearance. The highest confusion among species/genera found on different continents in TD2 was Coluber constrictor as Hierophis viridiflavus (7%). Both are slender, fast-moving, diurnal snakes, and the subspecies H. v. carbonarius (sometimes considered a full species; Mezzasalma et al., 2015) from Italy is all black, making it very similar in appearance to adult forms of C. c. constrictor and C. c. priapus in the eastern United States.

Species pairs that were commonly confused in both TD1 and TD2 were Coluber constrictor as Pantherophis alleghaniensis (9%; also confused 6% of the time in TD1), Opheodrys vernalis as Opheodrys aestivus (8%; also confused 23% of the time in TD1) and Nerodia fasciata as Nerodia erythrogaster (7%; also confused 11% of the time in TD1).

In our third test dataset (TD3, containing images for just 27 of the 45 classes hand-selected from other sources), F1 = 0.53 and log-loss = 1.19 when all 45 classes were considered (TD3a), and F1 = 0.73 and log-loss = 1.03 when only the 27 classes relevant to both datasets were considered (TD3b; Table 1). In TD3, species classification accuracy ranged from 100% for 10 species to 25% for Nerodia erythrogaster. Interclass confusion was 0 for 648/729 (89%) of class pairs.

Among species that are uncontroversially delimited, there were 23 pairs with confusion >10% in TD3 (Table 3), including two that were also commonly confused in both TD1 and TD2 (Opheodrys vernalis as Opheodrys aestivus; 33% in TD3, 8% in TD2, 23% in TD1; Nerodia fasciata as Nerodia erythrogaster; 25% in TD3, 7% in TD2, 11% in TD1) and one that was also commonly confused in TD1 (Storeria occipitomaculata as Storeria dekayi; 18% in TD3 vs 10% in TD1). The highest confusion among species in different genera was 25% for two species pairs (Heterodon platirhinos as Nerodia erythrogaster and Storeria occipitomaculata as Nerodia erythrogaster) and the highest confusion among species that do not occur in sympatry was Charina bottae as Haldea striatula (22%) (see Table 3 for detailed comparison).

When we asked human experts to identify the same images, they performed slightly better overall (species-level F1 = 0.76 for 45 classes, F1 = 0.79 for 27 classes; Table 1), with significant variation among species (Durso et al., 2021). Humans correctly labeled the 248 images in TD3 at the species level 53% of the time (N = 26,672), although this varied by species from 83% for Agkistrodon contortrix to 30% for Nerodia erythrogaster.

Thirteen species pairs were confused >10% of the time by humans, five of which also exceeded 10% in at least one of the three test datasets and all but one of which involved species that were also commonly confused by the algorithm (Table 3). Among frequently confused species from TD1 and TD2, humans only had the option to select “Black Ratsnake complex”, so there were no opportunities for confusion among the three putative species (P. obsoletus, P. alleghaniensis, and P. spiloides). Recent changes to taxonomy resulting in unfamiliar names may also hinder human experts’ ability to identify snakes (Carrasco et al., 2016), which could partially explain why e.g., Haldea striatula (more widely known as Virginia striatula in recent decades; Powell et al., 1994; McVay and Carstens, 2013) was correctly identified to the species level by humans only 31% of the time, or why Lampropeltis californiae was identified as Lampropeltis getula (the species from which it was split; Pyron and Burbrink, 2009) 23% of the time.

Some of the images in this dataset were purposefully chosen to be as difficult to identify as possible—e.g., low resolution, out of focus, with the snake filling only a small part of the frame, and/or obscured by vegetation. Stallkamp et al. (2012) also found that humans were superior to computer vision at identifying images with visual artifacts (e.g., traffic signs with graffiti or a lot of glare).

Because TD1 is a withheld subset of the training data and TD2 was collected using very similar methods, their similarity to the training data is very high. In contrast, TD3 images were selected to be representative of realistic difficult cases. We suggest that this explains the drop in performance for TD3, because image datasets often have built-in biases that are difficult to pinpoint and tests of cross-dataset generalization are rare (Dollár et al., 2009) but generally show that algorithms perform better on their “native” test sets (here, TD1 and TD2) than on entirely novel testing datasets (here, TD3) (Torralba and Efros, 2011). Of the types of bias discussed by Torralba and Efros (2011), we took pains to eliminate or account for label bias (inconsistent category definitions across datasets) by using the same snake taxonomy consistently and paying particular attention to cases where instability in species definitions might sabotage identification predictions. Our dataset is probably most susceptible to capture bias (photographers tending to take pictures of objects in similar ways) and perhaps to selection bias (source datasets that prefer particular kinds of images). Future studies should address the trade-offs between effort and information gain that may result from attempting to acquire images taken from particular angles or of particular anatomical features (Rzanny et al., 2017; Rzanny et al., 2019).

Species that were consistently identified with high accuracy across the three test datasets include:

• Boa constrictor (Boa Constrictor) had the highest average F1 (Figure 6), the fourth highest average recall (Supplementary Figure S3) and the third highest average precision (Supplementary Figure S4), as well as among the lowest average false positive (6.3 ± 2.7%) and average false negative (4.0 ± 2.3%) rates, across the three test datasets. These large, iconic snakes are easily recognized and there are few similar species with which they could be confused, especially within their range. Potential for confusion with large Python and Eunectes (anaconda) species exists. Several recent studies have suggested that Boa constrictor is likely to be a species complex containing as many as nine species (see Reynolds and Henderson, 2018 for a summary).

Species groups that were consistently confused across the three test datasets include:

• The two greensnake species in the genus Opheodrys: O. aestivus (Rough Greensnake) and O. vernalis (Smooth Greensnake; formerly placed in the genus Liochlorophis). Average F1 across test datasets and species was 81.4 ± 6.9% (Figure 6). These sister taxa are similar in their bright green color, slender bodies, small size, and gestalt; they differ anatomically mainly in that O. aestivus has keeled dorsal scales in 17 rows, whereas O. vernalis has smooth dorsal scales in 15 rows, subtle characteristics that may not be discernible in photos that are not taken at close range. Habitat and geographic range are also useful in distinguishing them (Ernst and Ernst, 2003). Although telling Opheodrys species apart is not likely to ever be clinically or epidemiologically useful, there are many other solid bright green snakes lacking easy-to-distinguish features in Africa, including non-venomous Philothamnus and Hapsidophrys as well as some highly dangerous Dispholidus (boomslangs) and Dendroaspis (mambas); see e.g., Broadley and Fitzsimons (1983) and Chippaux and Jackson (2019).

Of the 44 species pairs that exceeded 10% confusion on any test dataset (including human IDs), only four had no overlap in geographic range (two others had very little geographical overlap). Incorporating information on the geography has the potential to better discriminate these species pairs (Wittich et al., 2018), but in the most biodiverse regions of the world, >125 species of snakes may be sympatric within the same 50 × 50 km² area (Roll et al., 2017).

Most commonly confused species pairs involved either two non-venomous or two medically important venomous species (hereafter “venomous”). Of the 44 species pairs that exceeded 10% confusion on any test dataset (including human IDs), 31 involved two non-venomous species and five involved two venomous species (Table 3). Only three of the 44 pairs were “false positives” (a non-venomous species identified as a venomous one) and only five were “false negatives” (a venomous species identified as a non-venomous one). Bites from all medically important venomous snakes in our dataset (all North American pit vipers) would be treated using the same antivenom (Dart et al., 2001; Gerardo et al., 2017) but this would not be as straightforward in other regions of the world (e.g., Bitis vs. Echis in sub-Saharan Africa) and indeed is becoming more complex in North America as new products enter the market (Bush et al., 2015; Cocchio et al., 2020).

Finally, we found that among the 163 confused species pairs in TD3, the algorithm suggested the species with more images in the training dataset as the identity of the species with fewer images far more often (78% of pairs) than the other way around. For example, Opheodrys vernalis (1,471 training images) was misidentified as O. aestivus (3,312 training images) 33% of the time in TD3, but O. aestivus was never misidentified as O. vernalis in TD3, even though humans confused these two species nearly equally often in both directions (12 and 10%). Humans were marginally more likely to suggest the species from a given pair with more images even when it was incorrect (65 of 118 confused pairs; 55%), but the effect was much more pronounced for AI than for humans (Figure 7). For the algorithm, this is probably due to the imbalance in the amount of training data among classes; for humans, we suggest that it may be caused by the same processes that lead a species to be represented by fewer images in the training data: smaller geographic ranges or lower probabilities of encounter, leading to less familiarity with the rarer species.

Long-tailed distributions are commonplace and have been called “the enemy of machine learning” (Bengio, 2015). As we enlarge our dataset to include more of the snake species with few images per species (i.e., those in the “long-tailed” part of the distribution), we must pay particular attention to rare species. Previous studies on faces (Zhang et al., 2017) and scene parsing (Yang et al., 2014) suggest approaches including rare class enrichment (Yang et al., 2014) and clustering objects into visually similar hierarchical groups (in our case, snake genera or families) (Ouyang et al., 2016).