Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • spatialite
2 results
Blame Permalink
ExportToWorldShape.py 5.10 KiB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123 






#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 14 13:21:12 2021

Export the stationlist to the geopandas world shape file

@author: geopeter
"""
import numpy as np
import geopandas as gpd


def buildAverageTimeseries(stationList, fromYear, toYear, name):
    meanAverage = []
    for stationID, station in stationList.iterrows():
        temps = []
        for i in range(fromYear,toYear):
            if not np.isnan(station[i]):
                temps.append(station[i])
        
        if len(temps) > 5:
            meanAverage.append(np.mean(temps))
        else:
            meanAverage.append(np.NaN)
    
    stationList[name] = np.round(meanAverage,1)


def cleanAverageTimeseries(stationList):
    # determine gauges that includes both timeseries. If not delete them.
    
    for stationID, station in stationList.iterrows():
        if np.isnan(station['m1961T1990']) or np.isnan(station['m1991T2018']):
            #station['m1961T1990'] = None
            #station['m2010T2018'] = None
            stationList.at[stationID, "m1961T1990"] = None
            stationList.at[stationID, "m1991T2018"] = None

def renameCountryNames(stationList):
    
    # rename countries to merge with geopandas world shape file
    stationList.loc[stationList['country']=="Korea, Dem. People's Rep.", 'country'] = 'South Korea'
    stationList.loc[stationList['country']=="Slovakia (Slovak. Rep.)", 'country'] = 'Slovakia'
    stationList.loc[stationList['country']=="Slowenia", 'country'] = 'Slovenia'
    stationList.loc[stationList['country']=="Russian Federation", 'country'] = 'Russia'
    stationList.loc[stationList['country']=="Bosnia and Herzegowina", 'country'] = 'Bosnia and Herz.'
    stationList.loc[stationList['country']=="Slovakia (Slovak. Rep.)", 'country'] = 'Slovakia'
    stationList.loc[stationList['country']=="Croatia/Hrvatska", 'country'] = 'Croatia'
    stationList.loc[stationList['country']=="Moldova, Rep. Of", 'country'] = 'Moldova'
    stationList.loc[stationList['country']=="United Kingdom of Great Britain and N.-Ireland ", 'country'] = 'United Kingdom'
    stationList.loc[stationList['country']=="Czech Republic", 'country'] = 'Czechia'
    stationList.loc[stationList['country']=="Somalia", 'country'] = 'Somaliland'
    stationList.loc[stationList['country']=="Iran (Islamic Rep. of)", 'country'] = 'Iran'
    stationList.loc[stationList['country']=="Mauretania", 'country'] = 'Mauritania'
    stationList.loc[stationList['country']=="Central African Republic", 'country'] = 'Central African Rep.'
    stationList.loc[stationList['country']=="South Sudan", 'country'] = 'S. Sudan'
    stationList.loc[stationList['country']=="Dem. Republic of the Congo", 'country'] = 'Dem. Rep. Congo'
    stationList.loc[stationList['country']=="Mauretania", 'country'] = 'Somalia'
    stationList.loc[stationList['country']=="Syrian Arab Rep.", 'country'] = 'Syria'
    stationList.loc[stationList['country']=="Australien, SW-Pazifik", 'country'] = 'Australia'
    stationList.loc[stationList['country']=="Western-Sahara",'country'] = "W. Sahara"
    
    return stationList

def convertStationListToGPD(stationList):
    
    print("convert stationlist to GeoPandas")
    
    del stationList["file"]
    
    stationList = stationList[stationList.lat != ""]
    stationList = stationList[stationList.lon != ""]
    
    stationList = renameCountryNames(stationList)
    
    stationList["lat"] = stationList["lat"].astype(str).astype(float)
    stationList["lon"] = stationList["lon"].astype(str).astype(float)
    
    # add some average timeseries
    buildAverageTimeseries(stationList, 1961,1990,"m1961T1990")
    buildAverageTimeseries(stationList, 1991,2018,"m1991T2018")
    cleanAverageTimeseries(stationList)
    
    # convert stationlist to geodataframe
    stationGPD = gpd.GeoDataFrame(stationList, geometry=gpd.points_from_xy(stationList.lat, stationList.lon)).reset_index()
    # del stationGPD["lat"]
    # del stationGPD["lon"]
    stationGPD = stationGPD.set_crs('epsg:4326')
    
    stationGPD.columns = stationGPD.columns.astype(str)
    stationGPD = stationGPD.sort_index(axis=1, ascending=False)
    #stationGPD.to_file("stationList.shp", "ESRI Shapefile")
    
    return stationGPD
    
    
def exportStationListToWorldShape(stationGPD):
    
    print("export stationlist to world shape file")
    
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    countryMeanGPD = gpd.sjoin(world, stationGPD, how="inner", op='intersects')
    countryMeanGPD = countryMeanGPD.groupby("name", axis=0).mean().reset_index()
    
    countryMeanGPD["anomalie"] = countryMeanGPD["m1991T2018"] - countryMeanGPD["m1961T1990"]
    
    del countryMeanGPD["pop_est"]
    del countryMeanGPD["gdp_md_est"]
    
    for i in range(1873,1950):
        try:
            del countryMeanGPD[str(i)]
        except:
            # do nothing
            continue
    
    worldGauge = world.set_index("name").join(countryMeanGPD.set_index("name"))
    worldGauge.to_file("./output/countryAnnualTemperature.shp", "ESRI Shapefile")
    
def export(stationList):
    stationGPD = convertStationListToGPD(stationList)
    exportStationListToWorldShape(stationGPD)