472-final-project

This project is to cooperate with Agrocorp International LTD to make an infographic for them. When I start this project, I have considered which software is good and suitable for this map. The primal idea of mine is using R studio to make an online dynamic map for the company. According to my experience on assignment3, I find that R studio is good to process huge data and operate, but it also requires more knowledge for programming. My knowledge to R studio is very basic and limited, and I have realized the R will be difficult for me to make such global trade map. In addition, the dynamic map maybe better to show the process of global trade by ships, but it is not clear to see the market scale on the trade. Then I consider to use ArcGIS to make this map. The ArcGIS is good at data processing and analyzing, and map display, but it is not good-looking. As a result, I email to the company and listen to their suggestions and requirements for the map. According to the company’s interest and requirement, they want to have more focus on the top market and top products of export. Also, the company wants to use this map for their marketing meeting, so they require that map should be easy to read and understand, and the key information should be highlighted. Thus, I decide to make an infographic because it is a visual design that help to explain complicated data in a simple way. I have also considered using ArcGIS for the base map, and it is good to show the world. Then I will mainly use Adobe Illustrator for the infographic’s background and visual design. To do a good infographic, the data processing is primary important. First is the gathering data. After getting the raw data, next step is mining and filtering data. After collecting data, I start to do some calculations for the sum of agri-food by countries and by types. After the data processing, then I can start for the map design. The company requires an intuitionistic visual performance, so I decide to use a world map with Winkel tripel projection.To differentiate Canada with other countries, I use a larger font size to label Canada. And now, Canada is the most striking on the map, and can attract more attention. Next, I plot the other countries’ points at the location of their capital. Then, I connect each point to Canada to make lines. These lines will give people a sense of global trade, like the ship line. After finishing the world trade map, the background color is also important for the map design. From my idea, this infographic is talking about the agricultural products and overseas ship trade. Thus, I decide to choose green-blue gradient color to set up the background.After finishing the background and the world trade map, the last thing I have to do is transform data into beautiful graphics.For the font style, I use Kohinoor Bangla for all words and numbers because of its concision and clarity.

Link:https://boruizhang1992.carto.com/viz/238a8f9c-a60b-11e6-847f-0e98b61680bf/public_map

Deconstruction Critique “Chinese GDP equivalents” interactive map

Link: http://www.economist.com/content/chinese_equivalents

I choose the interactive web map from The Economist website. It compares Chinese provinces with other countries on GDP, GDP per person, population, and exports. This map shows the nearest equivalent countries to each province. First, the map designer acquired the original GDP and population data for each Chinese province and each country in the world. These data may be from WTO, China’s National Statistics Department and CEIC. Furthermore, to see the data clearly and conveniently, the designer processed these data to correct to one decimal. Then, he used the GDP to divide by population in order to get GDP per person. However, some data may be not credible. Some Chinese provinces may exaggerate their output: the sum of their reported GDPs is 10% higher than the national total. The designer also amended these data to let the overstatement is modest. After parsing and filtering data, he matched Chinese provinces to the nearest countries in order to make the “country equivalents” interactive map.

There are total four maps which focus on four different indexes. The first one is about GDP, the second one talks about GDP per person, the third one is focused on population, and the last one emphasizes exports. For each map, the designer used five colors to divide five levels from low to high. Then he labeled each relevant country on each Chinese province directly. When users click on each province, they can get detailed information about this province and applicable country on the same index. On the right side, there is a histogram distribution to show the ranking for each Chinese province. Different index brings different ranking, for example, even though Inner Mongolia sits in the middle of GDP ranking, it has a high GDP per person ranking. This histogram is better and clearer to show the regional inequality between each province.

The map designer compared the economy of the Chinese provinces to other world’s countries. His aim may be to highlight China’s economic power. The data for this map were collected in 2010 and this map was published in 2011. As we all know, during this period, many countries were experiencing a global economic crisis. On the contrary, China still kept the fast economic growth. Then lots of people including economists, government officials and this map designer, may have an optimistic expectation for future Chinese economy. Some of them even advocate China will replace the US to become the strongest economic country in the world. However, after 2014, the speed of Chinese economic development started to slow down. Many problems arose and were enlarged. For example, the gap of wealth, the structure of economic development sector, and the local governmental debt. Actually, these problems had occurred already in 2010, but the designer did not pay more attention to these problems. Thus, this GDP contrastive map cannot response the real economic development of China at that time, it only can show China’s huge economic quantity. In other words, this GDP per capital adjusted by PPP (purchasing power parity) cannot reflect China’s true economic power. Thus, this map’s representation does not link the reality.

To mend this map better in order to let the public know the real Chines economic power rationally. I think the designer should consider each Chinese province’s governmental debt and effective economic investment. GDP figures are useful but are often misleading. The figure that matters is the net “profits” an economy generates and not just the market value of what it sells. In addition, when taking Chinese province to compare to other countries, this metaphorical shows powerful Chinese economy and good development prospect. From my opinion, it is better and more rational to compare Chinese province to other countries’ states. Finally, for me, this kind of comparison is not meaningful, because after all, the average in china is still very low and there are still lots of people in poor. These attractive economic data cannot cover Chinese economic development problems. While if contrast the imbalance in regional development in China, it will be more significant.

# ------------------------------------------------------------------ #
# ---------------------- Install Libraries ------------------------- #
# ------------------------------------------------------------------ #
install.packages("GISTools")
install.packages("RJSONIO")
install.packages("rgdal")
install.packages("RCurl")
install.packages("curl")

# Unused Libraries:
# install.packages("ggmap")
# library(ggmap)
# ------------------------------------------------------------------ #
# ----------------------- Load Libararies -------------------------- #
# ------------------------------------------------------------------ #
library(GISTools)
library(RJSONIO)
library(rgdal)
library(RCurl)
library(curl)

Step1: This step is the beginning, and to make sure the R studio has installed and loaded the tools for my project.

# ------------------------------------------------------------------ #
# ---------------------------- Acquire ----------------------------- #
# ------------------------------------------------------------------ #
# access csv data from desktop
fname = ('/User/Borui/Desktop/crime_2004.csv')
# Read data as csv
data = crime_2004.csv(fname, header=T)
# inspect your data
print(head(crime_2004))

Step2: When the R studio is ready, then I import the data about Vancouver crime statistics in 2004 into R, and let it read from the local.

# ------------------------------------------------- #
# -------------- Parse: Geocoder ------------------- #
# ------------------------------------------------- #
# change intersection to 00's
crime_2004$HUNDRED_BLOCK = gsub("X","0",crime_2004$HUNDRED_BLOCK)
print(head(crime_2004$HUNDRED_BLOCK))

# Join the strings from each column together & add "Vancouver, BC":
crime_2004$full_address = paste(crime_2004$HUNDRED_BLOCK, 
 paste(crime_2004$NEIGHBOURHOOD,
 "Vancouver, BC",
 sep=", "),
 sep=" ")

# removing "Intersection " from the full_address entries
crime_2004$full_address = gsub("Intersection ", "", crime_2004$full_address)
print(head(crime_2004$full_address))

Step3: Then I notice some number ‘0’ show the symbol ‘X’ in the data. I change intersection to 00’s, and let ‘X’ to become ‘0’ under the HUNDRED_BLOCK column. Next, I created a new column called ‘full address’ to combine the HUNDRED_BLOCK with NEIGHBOURHOOD. Also I remove the intersection in order to make my data clearer and easier.

# BC geocoading
bc_geocode = function(search){
 if(!is.character(search)){stop("'search' must be a character string")}
 search = RCurl::curlEscape(search)
 base_url = "http://apps.gov.bc.ca/pub/geocoder/addresses.json?addressString="
 url_tail = "&locationDescriptor=any&maxResults=1&interpolation=adaptive&echo=true&setBack=0&outputSRS=4326&minScore=1&provinceCode=BC"
 final_url = paste0(base_url, search, url_tail)
 response = RCurl::getURL(final_url)
 response_parsed = RJSONIO::fromJSON(response)
 if(length(response_parsed$features[[1]]$geometry[[3]]) > 0){
 geocoords = list(lon = response_parsed$features[[1]]$geometry[[3]][1],
 lat = response_parsed$features[[1]]$geometry[[3]][2])
 }else{
 geocoords = NA
 }
 return(geocoords)
}

# Create an empty vector for latitude and longtitude coordinates
lat = c() 
lon = c()
# loop through the addresses
for(i in 1:length(crime_2004$full_address)){
 # store the address at index "i" as a character
 address = crime_2004$full_address[i]
 # append the latitude of the geocoded address to the lat vector
 lat = c(lat, bc_geocode(address)$lat)
 # append the longitude of the geocoded address to the lon vector
 lon = c(lon, bc_geocode(address)$lon)
 # at each iteration through the loop, print the coordinates
 print(paste("#", i, ", ", lat[i], lon[i], sep = ","))
}

# add the lat lon coordinates to the dataframe
crime_2004$lat = lat
crime_2004$lon = lon
# export csv file
ofile = "/Users/Borui/Desktop/crime.csv" 
write.csv(crime_2004, ofile)

Step4: This step is a huge and time-consuming task. I firstly define the bc geocoding function in order to run the bc geocoding API to get the latitude and longitude coordinates for each full address. This task costs me 4 hours. After a major computational task, my data now has specific latitude and longitude coordinates related to the crime scene. Then I output this data and named crime_csv to keep a copy.

# ------------------------------------------------- #
# --------------------- Mine ---------------------- #
# ------------------------------------------------- #
# --- Examine the unique cases --- #
unique(crime_2004$TYPE)
for (i in 1:length(unique(crime_2004$TYPE))){
 print(unique(crime_2004$TYPE)[i], max.levels=0)
}

#Divide classes to group type
#-----------------Mischief--------------#
Mischief = c('Mischief')
#-----------------Theft---------------#
Theft = c('Theft from Vehicle',
 'Theft of Vehicle',
 'Other Theft')
#----------------Break and Enter----------#
Break_and_Enter = c('Break and Enter Commercial',
 'Break and Enter Residential/Other')
#----------------Offence Against a Person------------#
Offence_Against_a_Person = c('Offence Against a Person')
#---------------Homicide--------------#
Homicide = c('Homicide')

# give class id numbers:
crime_2004$cid = 9999
for(i in 1:length(crime_2004$TYPE)){
 if(crime_2004$TYPE[i] %in% Mischief){
 crime_2004$cid[i] = 1 
 }else if(crime_2004$TYPE[i] %in% Theft){
 crime_2004$cid[i] = 2 
 }else if(crime_2004$TYPE[i] %in% Break_and_Enter){
 crime_2004$cid[i] = 3 
 }else if(crime_2004$TYPE[i] %in% Offence_Against_a_Person){
 crime_2004$cid[i] = 4 
 }else if(crime_2004$TYPE[i] %in% Homicide){
 crime_2004$cid[i] = 5 
 }else{
 crime_2004$cid[i] = 0 
 }
}

Step5: To find every crime types, I type a unique code. Then I integrate my data and divide to 5 different types crime. They are the mischief, the theft, the break and enter crime, offence against a person and the homicide. After the integration, I give the class ID number for each type crime.

# --- handle overlapping points --- #
# Set offset for points in same location:
crime_2004$lat_offset = crime_2004$lat
crime_2004$lon_offset = crime_2004$lon

# Run loop - if value overlaps, offset it by a random number
for(i in 1:length(crime_2004$lat)){
 if ( (crime_2004$lat_offset[i] %in% crime_2004$lat_offset) && (crime_2004$lon_offset[i] %in% crime_2004$lon_offset)){
 crime_2004$lat_offset[i] = crime_2004$lat_offset[i] + runif(1, 0.0001, 0.0005)
 crime_2004$lon_offset[i] = crime_2004$lon_offset[i] + runif(1, 0.0001, 0.0005)
 } 
}
# get a frequency distribution of the calls:
top_calls = data.frame(table(crime_2004$TYPE))
top_calls = top_calls[order(top_calls$Freq), ]
print(top_calls)

Step6: To avoid the same coordinates for the crime, I set the offset on the latitude and longitude by adding 0.0001 to 0.0005. Then I get two new columns for the lat_offset and lon_offset.

# ------------------------------------------------------------------ #
# ---------------------------- Filter ------------------------------ #
# ------------------------------------------------------------------ #
data_filter = subset(crime_2004, (lat <= 49.313162) & (lat >= 49.199554) & 
 (lon <= -123.019028) & (lon >= -123.271371) & is.na(lon) == FALSE )
# plot the data
plot(data_filter$lon, data_filter$lat)
ofile_filtered = '/Users/Borui/Desktop/crime_filtered.csv' 
write.csv(data_filter, ofile_filtered)

Step7: After the offset the latitude and longitude, I need to filter some data outside Vancouver’s boundary. Then I give a limit to the boundary coordinate to only focus on Vancouver. I plot each data to this limited area and export as crime_filtered.csv file.

#-------------Saving Dataset------------#
# --- Convert Data to Shapefile --- #
# store coordinates in dataframe
coords_crime_2004 = data.frame(data_filter$lon_offset, data_filter$lat_offset)
# create spatialPointsDataFrame
data_shp = SpatialPointsDataFrame(coords = coords_crime_2004, data = data_filter)
# set the projection to wgs84
projection_wgs84 = CRS("+proj=longlat +datum=WGS84")
proj4string(data_shp) = projection_wgs84
# set an output folder for our geojson files:
geofolder = "/Users/Borui/Desktop/GEOB_472/"
# Join the folderpath to each file name:
opoints_shp = paste(geofolder, "calls_", "1401", ".shp", sep = "")
print(opoints_shp)# print this out to see where your file will go & what it will be called
opoints_geojson = paste(geofolder, "calls_", "1401", ".geojson", sep = "") 
print(opoints_geojson) # print this out to see where your file will go & what it will be called
# write the file to a shp
writeOGR(data_shp, opoints_shp, layer = "data_shp", driver = "ESRI Shapefile",
 check_exists = FALSE)
# write file to geojson
writeOGR(data_shp, opoints_geojson, layer = "data_shp", driver = "GeoJSON",
 check_exists = FALSE)

Step8: On this step, with a filtered dataset, I can write my file out to a shapefile and a geojson. I create latitude and longitute coordinates of the crime in WGS84 projection shape file by using filtered data.

#---------------------Creating Hexagon Grid-------------------#
# --- aggregate to a grid --- #
# ref: http://www.inside-r.org/packages/cran/GISTools/docs/poly.counts
# set the file name - combine the shpfolder with the name of the grid
grid_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m.geojson'

# read in the hex grid setting the projection to utm10n
hexgrid = readOGR(grid_fn, 'OGRGeoJSON')

# transform the projection to wgs84 to match the point file and store
# it to a new variable (see variable: projection_wgs84)
hexgrid_wgs84 = spTransform(hexgrid, projection_wgs84)

#-----------------Assigning Values to Hexagons-----------------#
# Use the poly.counts() function to count the number of occurrences of calls per grid cell
grid_cnt = poly.counts(data_shp, hexgrid_wgs84)

# define the output names:
ohex_shp = paste(geofolder, "hexgrid_250m_", "1401", "_cnts",
 ".shp", sep = "")
print(ohex_shp)

ohex_geojson = paste(geofolder, "hexgrid_250m_","1401","_cnts",".geojson") 
print(ohex_geojson)

# write the file to a shp
writeOGR(check_exists = FALSE, hexgrid_wgs84, ohex_shp, 
 layer = "hexgrid_wgs84", driver = "ESRI Shapefile")

# write file to geojson
writeOGR(check_exists = FALSE, hexgrid_wgs84, ohex_geojson, 
 layer = "hexgrid_wgs84", driver = "GeoJSON")

Step9: I download a hexgonal grid at a resolution of 250m and reproject it from UTM zone 10 to WGS84. Then using the polu.counts function to count the number of points in each grid cell.  Also output them to shapefile and geojson.

#---------------------Version1---------------------------#
#--------------------------------------------------------#
# add the leaflet library to your script
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# subset out your data_filter:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

# initiate leaflet
m = leaflet()
# add openstreetmap tiles (default)
m = addTiles(m)
# you can now see that our maptiles are rendered
m

colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'),
 domain = data_filter$cid)
m = addCircleMarkers(m, 
 lng = df_Mischief$lon_offset, # we feed the longitude coordinates 
 lat = df_Mischief$lat_offset,
 popup = df_Mischief$Case_Type, 
 radius = 2, 
 stroke = FALSE, 
 fillOpacity = 0.75, 
 color = colorFactors(df_Mischief$cid),
 group = "1 - Mischief")
m = addCircleMarkers(m, 
 lng = df_Theft$lon_offset, lat=df_Theft$lat_offset, 
 popup = df_Theft$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Theft$cid),
 group = "2 - Theft")
m = addCircleMarkers(m, 
 lng = df_Break_and_Enter$lon_offset, lat=df_Break_and_Enter$lat_offset, 
 popup = df_Break_and_Enter$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Break_and_Enter$cid),
 group = "3 - Break_and_Enter")
m = addCircleMarkers(m, 
 lng = df_Offence_Against_a_Person$lon_offset, lat=df_Offence_Against_a_Person$lat_offset, 
 popup = df_Offence_Against_a_Person$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Offence_Against_a_Person$cid),
 group = "4 - Offence_Against_a_Person")
m = addCircleMarkers(m, 
 lng = df_Homicide$lon_offset, lat=df_Homicide$lat_offset, 
 popup = df_Homicide$Case_Type, 
 radius = 2,
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Homicide$cid),
 group = "5 - Homicide")

m = addTiles(m, group = "OSM (default)") 
m = addProviderTiles(m,"Stamen.Toner", group = "Toner")
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite")

m = addLayersControl(m,
 baseGroups = c("Toner Lite","Toner"),
 overlayGroups = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide")
)
# make the map
m

Step10: This is the version1 code.  It is working with point data. Firstly, I use the subset function to get store my data into 5 different classes. Then I initiate the leaflet map and insure different colors for the different classes. Next I type addCircleMarker function to add points on the map.Finally, I use addProviderTiles function to change the map’s background and use addLayerControl function to create a control bar to choose the different layers.

#-------------------------Version2----------------------------#
#-------------------------------------------------------------#
# add the leaflet library to your script
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# subset out my crime_filtered:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

data_filterlist = list(df_Mischief = subset(data_filter, cid == 1),
 df_Theft = subset(data_filter, cid == 2),
 df_Break_and_Enter = subset(data_filter, cid == 3),
 df_Offence_Against_a_Person = subset(data_filter, cid == 4),
 df_Homicide = subset(data_filter, cid == 5))

# Remember we also had these groups associated with each variable? Let's put them in a list too:
layerlist = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide",
 "0 - other")

# We can keep that same color variable:
colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'), domain=data_filter$cid)

# Now we have our loop - each time through the loop, it is adding our markers to the map object:
for (i in 1:length(data_filterlist)){
 m = addCircleMarkers(m, 
 lng=data_filterlist[[i]]$lon_offset,
 lat=data_filterlist[[i]]$lat_offset, 
 popup=data_filterlist[[i]]$Case_Type, 
 radius=2,
 stroke = FALSE, 
 fillOpacity = 0.75,
 color = colorFactors(data_filterlist[[i]]$cid),
 group = layerlist[i]
 )
}
m = addTiles(m, "Stamen.TonerLite", group = "Toner Lite") 
m = addLayersControl(m,
 overlayGroups = c("Toner Lite"),
 baseGroups = layerlist
)
m

Step11: This is the version2 code. It is optimizing version. The same process like the version1, I also use subset function to divide my data. Then I create a list to include these variables. Other fundamental color and map code keep the same like the version1.

#----------------------Version3-1--------------------------#
#----------------------------------------------------------#
library(rgdal)
library(GISTools)
# define the filepath
hex_1401_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m_1401_counts.geojson'
# read in the geojson
hex_1401 = readOGR(hex_1401_fn, "OGRGeoJSON")

m = leaflet()
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite")
# Create a continuous palette function
pal = colorNumeric(
 palette = "Greens",
 domain = hex_1401$data
)

# add the polygons to the map
m = addPolygons(m, 
 data = hex_1401,
 stroke = FALSE, 
 smoothFactor = 0.2, 
 fillOpacity = 1,
 color = ~pal(hex_1401$data),
 popup = paste("Number of calls: ", hex_1401$data, sep="")
)
m = addLegend(m, "bottomright", pal = pal, values = hex_1401$data,
 title = "Vancouver Crime Distribution in 2004",
 labFormat = labelFormat(prefix = " "),
 opacity = 0.75
)
m

Step12: This is the version3-1 code. It is working with polygons. To do this map, I need to read my hegrid with my data. Then I use colorNumeric function to choose green to show on the map. Finally, I add my polygon on the map.

#---------------------------Version3-2----------------------------#
#-----------------------------------------------------------------#
library(rgdal)
library(GISTools)
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# initiate leaflet map layer
m = leaflet()
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite") 

# --- hex grid --- #
# store the file name of the geojson hex grid
hex_1401_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m_1401_counts.geojson'

# read in the data
hex_1401 = readOGR(hex_1401_fn, "OGRGeoJSON")

# Create a continuous palette function
pal = colorNumeric(
 palette = "Greens",
 domain = hex_1401$data
)

# add hex grid
m = addPolygons(m, 
 data = hex_1401,
 stroke = FALSE, 
 smoothFactor = 0.2, 
 fillOpacity = 1,
 color = ~pal(hex_1401$data),
 popup= paste("Number of calls: ",hex_1401$data, sep=""),
 group = "hex"
)

# add legend
m = addLegend(m, "bottomright", pal = pal, values = hex_1401$data,
 title = "Vancouver Crime Distribution in 2004",
 labFormat = labelFormat(prefix = " "),
 opacity = 0.75
)

# --- points data --- #
# subset out my crime_filtered:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

data_filterlist = list(df_Mischief = subset(data_filter, cid == 1),
 df_Theft = subset(data_filter, cid == 2),
 df_Break_and_Enter = subset(data_filter, cid == 3),
 df_Offence_Against_a_Person = subset(data_filter, cid == 4),
 df_Homicide = subset(data_filter, cid == 5))


layerlist = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide",
 "0 - other")


colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'), domain=data_filter$cid)
for (i in 1:length(data_filterlist)){
 m = addCircleMarkers(m, 
 lng=data_filterlist[[i]]$lon_offset, lat=data_filterlist[[i]]$lat_offset, 
 popup=data_filterlist[[i]]$Case_Type, 
 radius=2,
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(data_filterlist[[i]]$cid),
 group = layerlist[i]
 )
 
}

m = addLayersControl(m,
 overlayGroups = c("Toner Lite", "hex"),
 baseGroups = layerlist
)
# show map
m

Step13: This is another vesion3-2 code, but it is synthesis. It includes both points and polygons on the map. For the point part, I copy the code from version2, then I also keep the code from version3 in order to combine point and polygon maps.

Comments: At the beginning, I have around 56228 records. After running the data through passing and geocoding and filtering, I have 44918 records left. I think these data are reliable without lost of error. When seeing the final map, it is clear that the map only focus on Vancouver with different types of crime. A shpfile is a vector data storage format for storing the location, shape and attributes of geographic features. Csv means common separated values, and it is a plain text format. It holds plain text as a series of values. Excel or xls file holds information about all the worksheets in a workbook, comprising both content and formatting. Geojson is a standard for representing the core geographic features that we are interested in. It does not contain everything that a shapfile can contain and it does not contain any styling information. In this project, I use csv file (crime_2004.csv) for the data source, and use shpfile (calls_1401.shp) and geojson (hexgrid_250m.geojson) to store the location in order to show the the distribution on the final map.

---------------Full Script---------------
#####################################################################
# Vancouver Crime Distribution in 2004: Data Processing Script
# Date:2016.10.31
# By: Rui Zhang (Jerry)
#####################################################################
# ------------------------------------------------------------------ #
# ---------------------- Install Libraries ------------------------- #
# ------------------------------------------------------------------ #
install.packages("GISTools")
install.packages("RJSONIO")
install.packages("rgdal")
install.packages("RCurl")
install.packages("curl")

# Unused Libraries:
# install.packages("ggmap")
# library(ggmap)
# ------------------------------------------------------------------ #
# ----------------------- Load Libararies -------------------------- #
# ------------------------------------------------------------------ #
library(GISTools)
library(RJSONIO)
library(rgdal)
library(RCurl)
library(curl)

# ------------------------------------------------------------------ #
# ---------------------------- Acquire ----------------------------- #
# ------------------------------------------------------------------ #
# access csv data from desktop
fname = ('/User/Borui/Desktop/crime_2004.csv')
# Read data as csv
data = crime_2004.csv(fname, header=T)
# inspect your data
print(head(crime_2004))

# ------------------------------------------------- #
# -------------- Parse: Geocoder ------------------- #
# ------------------------------------------------- #
# change intersection to 00's
crime_2004$HUNDRED_BLOCK = gsub("X","0",crime_2004$HUNDRED_BLOCK)
print(head(crime_2004$HUNDRED_BLOCK))

# Join the strings from each column together & add "Vancouver, BC":
crime_2004$full_address = paste(crime_2004$HUNDRED_BLOCK, 
 paste(crime_2004$NEIGHBOURHOOD,
 "Vancouver, BC",
 sep=", "),
 sep=" ")

# removing "Intersection " from the full_address entries
crime_2004$full_address = gsub("Intersection ", "", crime_2004$full_address)
print(head(crime_2004$full_address))

# BC geocoading
bc_geocode = function(search){
 if(!is.character(search)){stop("'search' must be a character string")}
 search = RCurl::curlEscape(search)
 base_url = "http://apps.gov.bc.ca/pub/geocoder/addresses.json?addressString="
 url_tail = "&locationDescriptor=any&maxResults=1&interpolation=adaptive&echo=true&setBack=0&outputSRS=4326&minScore=1&provinceCode=BC"
 final_url = paste0(base_url, search, url_tail)
 response = RCurl::getURL(final_url)
 response_parsed = RJSONIO::fromJSON(response)
 if(length(response_parsed$features[[1]]$geometry[[3]]) > 0){
 geocoords = list(lon = response_parsed$features[[1]]$geometry[[3]][1],
 lat = response_parsed$features[[1]]$geometry[[3]][2])
 }else{
 geocoords = NA
 }
 return(geocoords)
}

# Create an empty vector for latitude and longtitude coordinates
lat = c() 
lon = c()
# loop through the addresses
for(i in 1:length(crime_2004$full_address)){
 # store the address at index "i" as a character
 address = crime_2004$full_address[i]
 # append the latitude of the geocoded address to the lat vector
 lat = c(lat, bc_geocode(address)$lat)
 # append the longitude of the geocoded address to the lon vector
 lon = c(lon, bc_geocode(address)$lon)
 # at each iteration through the loop, print the coordinates
 print(paste("#", i, ", ", lat[i], lon[i], sep = ","))
}


# add the lat lon coordinates to the dataframe
crime_2004$lat = lat
crime_2004$lon = lon
# export csv file
ofile = "/Users/Borui/Desktop/crime.csv" 
write.csv(crime_2004, ofile)

# ------------------------------------------------- #
# --------------------- Mine ---------------------- #
# ------------------------------------------------- #
# --- Examine the unique cases --- #
unique(crime_2004$TYPE)
for (i in 1:length(unique(crime_2004$TYPE))){
 print(unique(crime_2004$TYPE)[i], max.levels=0)
}

#Divide classes to group type
#-----------------Mischief--------------#
Mischief = c('Mischief')
#-----------------Theft---------------#
Theft = c('Theft from Vehicle',
 'Theft of Vehicle',
 'Other Theft')
#----------------Break and Enter----------#
Break_and_Enter = c('Break and Enter Commercial',
 'Break and Enter Residential/Other')
#----------------Offence Against a Person------------#
Offence_Against_a_Person = c('Offence Against a Person')
#---------------Homicide--------------#
Homicide = c('Homicide')

# give class id numbers:
crime_2004$cid = 9999
for(i in 1:length(crime_2004$TYPE)){
 if(crime_2004$TYPE[i] %in% Mischief){
 crime_2004$cid[i] = 1 
 }else if(crime_2004$TYPE[i] %in% Theft){
 crime_2004$cid[i] = 2 
 }else if(crime_2004$TYPE[i] %in% Break_and_Enter){
 crime_2004$cid[i] = 3 
 }else if(crime_2004$TYPE[i] %in% Offence_Against_a_Person){
 crime_2004$cid[i] = 4 
 }else if(crime_2004$TYPE[i] %in% Homicide){
 crime_2004$cid[i] = 5 
 }else{
 crime_2004$cid[i] = 0 
 }
}

# --- handle overlapping points --- #
# Set offset for points in same location:
crime_2004$lat_offset = crime_2004$lat
crime_2004$lon_offset = crime_2004$lon

# Run loop - if value overlaps, offset it by a random number
for(i in 1:length(crime_2004$lat)){
 if ( (crime_2004$lat_offset[i] %in% crime_2004$lat_offset) && (crime_2004$lon_offset[i] %in% crime_2004$lon_offset)){
 crime_2004$lat_offset[i] = crime_2004$lat_offset[i] + runif(1, 0.0001, 0.0005)
 crime_2004$lon_offset[i] = crime_2004$lon_offset[i] + runif(1, 0.0001, 0.0005)
 } 
}
# get a frequency distribution of the calls:
top_calls = data.frame(table(crime_2004$TYPE))
top_calls = top_calls[order(top_calls$Freq), ]
print(top_calls)

# ------------------------------------------------------------------ #
# ---------------------------- Filter ------------------------------ #
# ------------------------------------------------------------------ #
data_filter = subset(crime_2004, (lat <= 49.313162) & (lat >= 49.199554) & 
 (lon <= -123.019028) & (lon >= -123.271371) & is.na(lon) == FALSE )
# plot the data
plot(data_filter$lon, data_filter$lat)
ofile_filtered = '/Users/Borui/Desktop/crime_filtered.csv' 
write.csv(data_filter, ofile_filtered)

#-------------Saving Dataset------------#
# --- Convert Data to Shapefile --- #
# store coordinates in dataframe
coords_crime_2004 = data.frame(data_filter$lon_offset, data_filter$lat_offset)

# create spatialPointsDataFrame
data_shp = SpatialPointsDataFrame(coords = coords_crime_2004, data = data_filter)

# set the projection to wgs84
projection_wgs84 = CRS("+proj=longlat +datum=WGS84")
proj4string(data_shp) = projection_wgs84

# set an output folder for our geojson files:
geofolder = "/Users/Borui/Desktop/GEOB_472/"

# Join the folderpath to each file name:
opoints_shp = paste(geofolder, "calls_", "1401", ".shp", sep = "")
print(opoints_shp)# print this out to see where your file will go & what it will be called

opoints_geojson = paste(geofolder, "calls_", "1401", ".geojson", sep = "") 
print(opoints_geojson) # print this out to see where your file will go & what it will be called

# write the file to a shp
writeOGR(data_shp, opoints_shp, layer = "data_shp", driver = "ESRI Shapefile",
 check_exists = FALSE)

# write file to geojson
writeOGR(data_shp, opoints_geojson, layer = "data_shp", driver = "GeoJSON",
 check_exists = FALSE)

#---------------------Creating Hexagon Grid-------------------#
# --- aggregate to a grid --- #
# ref: http://www.inside-r.org/packages/cran/GISTools/docs/poly.counts
# set the file name - combine the shpfolder with the name of the grid
grid_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m.geojson'

# read in the hex grid setting the projection to utm10n
hexgrid = readOGR(grid_fn, 'OGRGeoJSON')

# transform the projection to wgs84 to match the point file and store
# it to a new variable (see variable: projection_wgs84)
hexgrid_wgs84 = spTransform(hexgrid, projection_wgs84)

#-----------------Assigning Values to Hexagons-----------------#
# Use the poly.counts() function to count the number of occurrences of calls per grid cell
grid_cnt = poly.counts(data_shp, hexgrid_wgs84)

# define the output names:
ohex_shp = paste(geofolder, "hexgrid_250m_", "1401", "_cnts",
 ".shp", sep = "")
print(ohex_shp)

ohex_geojson = paste(geofolder, "hexgrid_250m_","1401","_cnts",".geojson") 
print(ohex_geojson)

# write the file to a shp
writeOGR(check_exists = FALSE, hexgrid_wgs84, ohex_shp, 
 layer = "hexgrid_wgs84", driver = "ESRI Shapefile")

# write file to geojson
writeOGR(check_exists = FALSE, hexgrid_wgs84, ohex_geojson, 
 layer = "hexgrid_wgs84", driver = "GeoJSON")

#---------------------Version1---------------------------#
#--------------------------------------------------------#
# add the leaflet library to your script
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# subset out your data_filter:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

# initiate leaflet
m = leaflet()
# add openstreetmap tiles (default)
m = addTiles(m)
# you can now see that our maptiles are rendered
m

colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'),
 domain = data_filter$cid)
m = addCircleMarkers(m, 
 lng = df_Mischief$lon_offset, # we feed the longitude coordinates 
 lat = df_Mischief$lat_offset,
 popup = df_Mischief$Case_Type, 
 radius = 2, 
 stroke = FALSE, 
 fillOpacity = 0.75, 
 color = colorFactors(df_Mischief$cid),
 group = "1 - Mischief")
m = addCircleMarkers(m, 
 lng = df_Theft$lon_offset, lat=df_Theft$lat_offset, 
 popup = df_Theft$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Theft$cid),
 group = "2 - Theft")
m = addCircleMarkers(m, 
 lng = df_Break_and_Enter$lon_offset, lat=df_Break_and_Enter$lat_offset, 
 popup = df_Break_and_Enter$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Break_and_Enter$cid),
 group = "3 - Break_and_Enter")
m = addCircleMarkers(m, 
 lng = df_Offence_Against_a_Person$lon_offset, lat=df_Offence_Against_a_Person$lat_offset, 
 popup = df_Offence_Against_a_Person$Case_Type, 
 radius = 2, 
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Offence_Against_a_Person$cid),
 group = "4 - Offence_Against_a_Person")
m = addCircleMarkers(m, 
 lng = df_Homicide$lon_offset, lat=df_Homicide$lat_offset, 
 popup = df_Homicide$Case_Type, 
 radius = 2,
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(df_Homicide$cid),
 group = "5 - Homicide")

m = addTiles(m, group = "OSM (default)") 
m = addProviderTiles(m,"Stamen.Toner", group = "Toner")
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite")

m = addLayersControl(m,
 baseGroups = c("Toner Lite","Toner"),
 overlayGroups = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide")
)

# make the map
m

#-------------------------Version2----------------------------#
#-------------------------------------------------------------#
# add the leaflet library to your script
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# subset out my crime_filtered:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

data_filterlist = list(df_Mischief = subset(data_filter, cid == 1),
 df_Theft = subset(data_filter, cid == 2),
 df_Break_and_Enter = subset(data_filter, cid == 3),
 df_Offence_Against_a_Person = subset(data_filter, cid == 4),
 df_Homicide = subset(data_filter, cid == 5))

# Remember we also had these groups associated with each variable? Let's put them in a list too:
layerlist = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide",
 "0 - other")

# We can keep that same color variable:
colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'), domain=data_filter$cid)

# Now we have our loop - each time through the loop, it is adding our markers to the map object:
for (i in 1:length(data_filterlist)){
 m = addCircleMarkers(m, 
 lng=data_filterlist[[i]]$lon_offset,
 lat=data_filterlist[[i]]$lat_offset, 
 popup=data_filterlist[[i]]$Case_Type, 
 radius=2,
 stroke = FALSE, 
 fillOpacity = 0.75,
 color = colorFactors(data_filterlist[[i]]$cid),
 group = layerlist[i]
 )
}
m = addTiles(m, "Stamen.TonerLite", group = "Toner Lite") 
m = addLayersControl(m,
 overlayGroups = c("Toner Lite"),
 baseGroups = layerlist
)
m

#----------------------Version3-1--------------------------#
#----------------------------------------------------------#
library(rgdal)
library(GISTools)
# define the filepath
hex_1401_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m_1401_counts.geojson'
# read in the geojson
hex_1401 = readOGR(hex_1401_fn, "OGRGeoJSON")

m = leaflet()
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite")
# Create a continuous palette function
pal = colorNumeric(
 palette = "Greens",
 domain = hex_1401$data
)

# add the polygons to the map
m = addPolygons(m, 
 data = hex_1401,
 stroke = FALSE, 
 smoothFactor = 0.2, 
 fillOpacity = 1,
 color = ~pal(hex_1401$data),
 popup = paste("Number of calls: ", hex_1401$data, sep="")
)
m = addLegend(m, "bottomright", pal = pal, values = hex_1401$data,
 title = "Vancouver Crime Distribution in 2004",
 labFormat = labelFormat(prefix = " "),
 opacity = 0.75
)
m

#---------------------------Version3-2----------------------------#
#-----------------------------------------------------------------#
library(rgdal)
library(GISTools)
library(leaflet)
filename = "C:/Users/Administrator/Desktop/crime_filtered.csv"
data_filter = read.csv(curl('https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/201401CaseLocationsDetails-geo-filtered.csv'), header=TRUE)
# initiate leaflet map layer
m = leaflet()
m = addProviderTiles(m, "Stamen.TonerLite", group = "Toner Lite") 

# --- hex grid --- #
# store the file name of the geojson hex grid
hex_1401_fn = 'https://raw.githubusercontent.com/joeyklee/aloha-r/master/data/calls_2014/geo/hgrid_250m_1401_counts.geojson'

# read in the data
hex_1401 = readOGR(hex_1401_fn, "OGRGeoJSON")

# Create a continuous palette function
pal = colorNumeric(
 palette = "Greens",
 domain = hex_1401$data
)

# add hex grid
m = addPolygons(m, 
 data = hex_1401,
 stroke = FALSE, 
 smoothFactor = 0.2, 
 fillOpacity = 1,
 color = ~pal(hex_1401$data),
 popup= paste("Number of calls: ",hex_1401$data, sep=""),
 group = "hex"
)

# add legend
m = addLegend(m, "bottomright", pal = pal, values = hex_1401$data,
 title = "Vancouver Crime Distribution in 2004",
 labFormat = labelFormat(prefix = " "),
 opacity = 0.75
)

# --- points data --- #
# subset out my crime_filtered:
df_Mischief = subset(data_filter, cid == 1)
df_Theft = subset(data_filter, cid == 2)
df_Break_and_Enter = subset(data_filter, cid == 3)
df_Offence_Against_a_Person = subset(data_filter, cid == 4)
df_Homicide = subset(data_filter, cid == 5)

data_filterlist = list(df_Mischief = subset(data_filter, cid == 1),
 df_Theft = subset(data_filter, cid == 2),
 df_Break_and_Enter = subset(data_filter, cid == 3),
 df_Offence_Against_a_Person = subset(data_filter, cid == 4),
 df_Homicide = subset(data_filter, cid == 5))

layerlist = c("1 - Mischief", "2 - Theft",
 "3 - Break_and_Enter","4 - Offence_Against_a_Person", "5 - Homicide",
 "0 - other")

colorFactors = colorFactor(c('red', 'orange', 'purple', 'blue', 'pink', 'brown'), domain=data_filter$cid)
for (i in 1:length(data_filterlist)){
 m = addCircleMarkers(m, 
 lng=data_filterlist[[i]]$lon_offset, lat=data_filterlist[[i]]$lat_offset, 
 popup=data_filterlist[[i]]$Case_Type, 
 radius=2,
 stroke = FALSE, fillOpacity = 0.75,
 color = colorFactors(data_filterlist[[i]]$cid),
 group = layerlist[i]
 )
 }
m = addLayersControl(m,
 overlayGroups = c("Toner Lite", "hex"),
 baseGroups = layerlist
)
# show map
m

Visualization: 

When making this map, I think the cartographic constraint is the final map, which combines the point with the polygon. When seeing this final map, there are so many points with deep colour, and the green polygons are covered by these points. As a result, it is unclear for people to see the frequency distribution on the map through the green polygons.

Insight:

According to the map, I can see downtown has the highest crime frequency no matter which types of crime. In addition, the theft and break and enter crime both are the most common rascality in Vancouver.  Downtown has the largest population density and lots of homeless people. It is make sense that these floating population bring more criminal activities. Furthermore, when focusing on downtown, I also notice that the east side has more criminal activities than the west side. Related to the real world, the east side is actually bad than other places, especially in the main street of downtown. Here is filled with lots of drug activities, erotic trades and drunken riots. There is no doubt that these bad activities will bring more criminal activities.

In addition to downtown, the distribution of theft and breaking crime in Vancouver area . When I see these criminal activities, they are more common on the busy street like the boardway and kingsway. Generally, the east side of Vancouver has more criminal activities than the west side of Vancouver. This is maybe caused by the population class. Due to the low rent and house price, there are more low educated and jobless people in the east side. These people maybe steal valuable things to sell. As a result, more theft and breaking crime happen here.

Reflection:

Now, I have experienced ArcGIS, Adobe Illustrator and R studio for the map making. My favourite software is the AI. It is good to design the infographic with a nice picture, but it is hard to deal with the enormous data. To process the huge data, R studio is better to operate. However, when comparing it to other softwares, it requires more base knowledge for programming. Thus, it has a high level difficulty for most people to design. The ArcGIS is good at data processing and map display. However, when comparing its map to the AI, it is not good-looking; when comparing its data processing to the R studio, it is not fast and accurate.

Comments on the front page of the UN Refugee Agency website.

From my point, the front page of the UN Refugee Agency is pretty good. First, the color choice is good and clear. The designer used light pink to represent the refugees, dark blue to show the asylum-seekers and light green to represent the internally displaced persons. The designer also used a light grey to show the world map. The combination of these colors not only makes a good contrast, but also increase readers’ interests on this map. Second, the user experience is also good. When people move the mouse over the map, they can notice the specific number occurring on the country. Then they click on the country, the chosen one will be highlight and other countries will become dim. This is more convenient for readers to choose countries and get the data. In addition, this map collects the refugee’s data from 1951 to 2015, and a histogram is created for the statistics. When readers click on other years, the corresponding map will change. This is easier for readers to see the differences during these years. However, the most obvious defect is the circle size. Even though the circle is enough big and clear for some countries, like Syrian, Iraq and Colombia, it is very small or invisible for some South America countries, like Bolivia, Argentina and Peru. Readers cannot see the circle occurring on these countries, but when they click these countries, the number of refugees is still existent. Moreover, the legend and explanations are clear and detailed for the readers. The designer made a particular interpretation for the representative colors. As a result, to be an infographic, this page contains good and visible map, detailed text information and clear graphs. On the whole, this front page is attractive and beautiful.

Original-data

Data-for-mine

Inforgraphic

Welcome to UBC Blogs. This is your first post. Edit or delete it, then start blogging!

size(600, 600);
background(186, 164, 131);

noStroke();
fill(196, 181, 161);
ellipse(300, 300, 500, 500);

noStroke();
fill(0, 0, 0);
quad(158.15, 305.7, 119.17, 369.68, 147.12, 347.62, 185.36, 284.37);
quad(236.67, 347.17, 168.43, 459.5, 216.67, 421.5, 285.33, 307.5);

noStroke();
fill(162, 153, 130);
quad(185.36, 284.37, 147.12, 347.62, 152.26, 351.25, 189.75, 289);
quad(285.33, 307.5, 216.67, 421.5, 221.67, 427.83, 291.33, 314.17);

stroke(188,106,84);
strokeWeight(5);
line(127.67, 414.5, 168, 463.83);
line(282.67, 284.37, 127.67, 414.5);

stroke(42,52,82);
strokeWeight(5);
line(216.67, 423.17, 170.33, 461.67);
line(156.25, 348.75, 216.67, 423.17);
line(264,259,156.25,348.75);

stroke(0);
strokeWeight(4);
noFill();
ellipse(338.86,261.06,10,10);
ellipse(339.53,259.17,22,22);

stroke(0);
strokeWeight(3.5);
noFill();
ellipse(341.61,256.21,38,38);

stroke(0);
strokeWeight(3);
noFill();
ellipse(343.08,253.28,52,52);

stroke(0);
strokeWeight(2.5);
noFill();
ellipse(345.27,250.9,64,64);
ellipse(346.83,247.83,76,76);
ellipse(348.67,244.67,88,88);
ellipse(350.59,241.34,100,100);
ellipse(352.5,238,112,112);
ellipse(353.83,234.34,124,124);
ellipse(356.33,231.75,136,136);
ellipse(357.83,229,148,148);
ellipse(359.67,226.17,160,160);
ellipse(362,222.83,172,172);
ellipse(363.67,219.5,184,184);
ellipse(365.33,216.67,196,196);
ellipse(366.83,214.5,208,208);

—————————————–

Remix Image:

void setup(){
size(600, 600);
background(67, 68, 68);

noStroke();
fill(165,165,165);
ellipse(303, 311, 524, 524);

stroke(130,130,130);
strokeWeight(8);
noFill();
ellipse(303,311,480,480);

stroke(196,196,196);
strokeWeight(6);
noFill();
ellipse(303,311,430,430);

stroke(107,110,167);
strokeWeight(4);
fill(107,110,167);
ellipse(303,311,10,10);

stroke(107,110,167);
strokeWeight(4);
noFill();
ellipse(303,311,22,22);

stroke(107,110,167);
strokeWeight(3.5);
noFill();
ellipse(303,311,34,34);

stroke(107,110,167);
strokeWeight(3);
noFill();
ellipse(303,311,46,46);

stroke(152,155,201);
strokeWeight(2.5);
noFill();
ellipse(303,311,58,58);
ellipse(303,311,70,70);
ellipse(303,311,82,82);
ellipse(303,311,94,94);
ellipse(303,311,106,106);
ellipse(303,311,118,118);

stroke(172,172,210);
strokeWeight(2);
ellipse(303,311,130,130);
ellipse(303,311,142,142);
ellipse(303,311,154,154);
ellipse(303,311,166,166);
ellipse(303,311,178,178);
ellipse(303,311,190,190);
ellipse(303,311,202,202);

stroke(188,190,222);
strokeWeight(1.5);
ellipse(303,311,214,214);
ellipse(303,311,226,226);
ellipse(303,311,238,238);
ellipse(303,311,250,250);
ellipse(303,311,262,262);
ellipse(303,311,274,274);
ellipse(303,311,286,286);
ellipse(303,311,298,298);

stroke(195,195,230);
strokeWeight(0.6);
ellipse(303,311,310,310);
ellipse(303,311,322,322);
ellipse(303,311,334,334);
ellipse(303,311,346,346);
ellipse(303,311,358,358);
ellipse(303,311,370,370);

stroke(231,29,24);
strokeWeight(2);
line(303,49,41,311);

stroke(250,238,0);
strokeWeight(2);
line(41,311,303,575);

stroke(34,172,56);
strokeWeight(2);
line(565,311,303,575);

stroke(46,167,224);
strokeWeight(2);
line(303,49,565,311);

}
void draw(){
  if(mousePressed){
    background(255);
  }
  ellipse(mouseX,mouseY,10,10);
}