Mapping population using tmap orders incorrectly

Question

I'm trying to map population in California counties for 2019. Using class() I saw the data is stored as characters so I changed it to numeric to get accurate results for summary(). When mapping, the values are being separated in an incorrect order (see legend). Not sure how else to solve this!

camap <- readOGR("CA_Counties_Pop2019.shp")
camap <- spTransform(camap, CRS("+init=epsg:4326"))

pop2019 <- as.numeric(as.character(camap$POP_2019))
summary(pop2019)

camap$poprange[pop2019 >= 1000 & pop2019 < 100000] <- "1,000-100,000"
camap$poprange[pop2019 >= 100000 & pop2019  < 500000] <- "100,000-500,000"
camap$poprange[pop2019 >= 500000 & pop2019  < 1000000] <- "500,000 - 1,000,000"
camap$poprange[pop2019 >= 1000000 & pop2019  < 5000000] <- "1,000,000 - 5,000,000"
camap$poprange[pop2019 >= 5000000 & pop2019  < 11000000] <- "5,000,000 - 11,000,000"

tm_shape(camap) + tm_fill("poprange", palette="Reds")`

• The map: