0

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-

  1. Lightning - 24 hrs period or annually
  2. Region ID
  3. 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)
2
  • 2
    Your code is a Jupyter code format ( *.ipynb file) and not a Python code. Use jupyter nbconvert --to script your_nodebook.ipynbto convert it to a Python code please – gene Feb 19 at 10:36
  • What's this "GHRC" data? You've not given us a link to it. – Spacedman Feb 19 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.