My code which I have effectively takes in as input the lightning data (in its case- http://wwlln.net/ but in my case, it's GHRC data), the global shapefiles, and population density data. How would I etch the GHRC data instead of wwlln? If I change it, what other parameters I also need to change?
**My three output would be-
- Lightning - 24 hrs period or annually
- Region ID
- Date**
I have attached the code below:
import geopandas as gpd
from shapely.geometry import Point
from rasterstats import point_query
import fiona
import pandas as pd
import gzip
import os
import re
from datetime import datetime
'''
:param path: Root folder for the data
:param outputfolder: subfolder of 'WLLN_Lightning/output/' where the data is saved (need to be manually created)
:param admin_file1: shapefile with the boundaries of the GWP subnational regions, will be appended with borders_2
:param admin_file2: shapefile with the boundaries of the GWP subnational regions, will be appended with borders_2
:param admin_id: variable used as admin_id for borders_1 + borders_2
:param pop_density: the raster from which population should be measured
:param ymin, ymax: year range
overwrite = False will skip days already computed
'''
path = '/Users/nikitamelnikov/Dropbox/WWLLN_Lightning/'
outputfolder = '3G_gov_weighted'
admin_file1 = 'admin/3G_gov/level1_map.shp'
admin_file2 = 'admin/3G_gov/level2_map.shp'
admin_id = 'regionid_m'
pop_density = 'admin/3G_gov/population_density.tif'
ymin, ymax = 2005, 2017
overwrite = False
os.chdir(path)
# GLOBALS (should not be modified)
##################################
# WWLLN
'''
Parameters for WWLLN data
corrupted_files are two files I could not read correctly on my computer,
it is manually fixed later in the script
'''
crs = {'init': 'epsg:4326'}
colnames = ['Date', 'Timestamp', 'Latitude', 'Longitude', 'Residuals', 'Nstations']
corrupted_files = ['A20120627.loc.gz', 'A20170713.loc.gz']
# ADMIN
'''
bounding box of the administrative limits,
the computation is made only for lightnings within that box
'''
borders_1 = gpd.read_file(path + admin_file1)
borders_2 = gpd.read_file(path + admin_file2)
borders = borders_1.append(borders_2, sort=False)
left, bottom, right, top = borders.total_bounds
# SCRIPT
########
for i, archive_name in enumerate(os.listdir(path+'rawdata/Afiles/')):
#for archive_name in ["A20080109.loc.gz"]:
# META
'''
check file name, create output name, and extract year
'''
assert archive_name[-7:] == '.loc.gz'
loc = path+'output/{}/geocoded_daily/{}.csv'
output_name = loc.format(outputfolder, archive_name[:-7])
year = int(archive_name[1:5])
# FILTER
'''
skip irrelevant files
'''
if os.path.isfile(output_name) and not overwrite:
continue
if year < ymin or year > ymax:
continue
print(archive_name) # just to see where we are in the compuation
# start = datetime.now() # just to see how long things take
# READ .GZ
with gzip.open(path+'rawdata/Afiles/' + archive_name , 'rb') as archive:
data = archive.read().decode('ascii')
data = [row.split(',') for row in data.split('\n')][:-1]
# print('Reading .gz:', datetime.now() - start)
# start = datetime.now()
# HANDLE CORRUPTED FILE (1)
if archive_name in corrupted_files:
if archive_name == 'A20120627.loc.gz':
data[286493] = ['2012/06/27', '12:44:58.552584', '-10.0520', ' 156.0042', '-99', '-99']
data[345053] = ['2012/06/27', '15:25:21.352359', ' -1.9076', ' -63.6534', '-99', '-99']
data += [['2012/06/27', '12:51:57.077721', '44.0976', '41.7293', '9.4', '5']]
data += [['2012/06/27', '12:51:57.077721', '44.0976', '41.7293', '9.4', '5']]
if archive_name == 'A20170713.loc.gz':
data[410561] = ['2017/07/13', '16:04:16.531394', ' 28.7434', ' -91', '-99', '-99']
data += [['2017/07/13', '16:10:00.000016', '5.3827', '100.2102', '12.1', '5']]
# TO PANDAS
data = pd.DataFrame(data, columns=colnames)
# HANDLE CORRUPTED FILE (2)
'''
just to remember it later if needed
'''
if archive_name in corrupted_files:
data['corrupted_file'] = 1
# FORMAT DATA
data.Longitude = data.Longitude.map(float)
data.Latitude = data.Latitude.map(float)
# print('Format:', datetime.now() - start)
# start = datetime.now()
# CLIP DATA (~)
data = data[data.Longitude >= left]
data = data[data.Longitude <= right]
data = data[data.Latitude >= bottom]
data = data[data.Latitude <= top]
# print('Clip:', datetime.now() - start)
# start = datetime.now()
if data.shape[0] > 0:
'''
only if some lightnings fall into the box that day
'''
# CONVERT TO GEODATAFRAME
geometry = [Point(x, y) for x, y in zip(data.Longitude, data.Latitude)]
lightnings = gpd.GeoDataFrame(data, crs=crs, geometry=geometry)
# print('Convert to geo:', datetime.now() - start)
# start = datetime.now()
# SPATIAL JOIN
out = gpd.sjoin(lightnings, borders, op='within', how='inner')
# print('Spatial join:', datetime.now() - start)
# start = datetime.now()
if out.shape[0] > 0:
'''
only if some lightnings fall into the admin limits that day
'''
# POP DENSITY
out.to_file(path+'tmp/tmp.shp')
stats = point_query(path+'tmp/tmp.shp', pop_density)
out['pop'] = stats
# print('Pop density:', datetime.now() - start)
# start = datetime.now()
# EXPORT (uncomment first line for full data)
#out[colnames + [admin_id, 'pop']].to_csv(output_name)
df = out.groupby(['regionid_m']).sum()
df.to_csv(output_name)
# print('Export:', datetime.now() - start)
jupyter nbconvert --to script your_nodebook.ipynb
to convert it to a Python code please – gene Feb 19 at 10:36