Bài 15: Phân tích Dữ liệu vector với GeoPandas

# Import các thư viện chính
import geopandas as gpd           # Thư viện chính cho phân tích địa không gian
import pandas as pd                # Xử lý dữ liệu bảng
from shapely.geometry import Point, Polygon
import os 
import warnings                  # Quản lý cảnh báo
warnings.filterwarnings('ignore')  # Ẩn cảnh báo không cần thiết
outpath = r'G:\My Drive\python\geocourse\data\vector'

geo_series = gpd.GeoSeries([
    Point(105.85, 21.02),  # Hà Nội
    Point(106.63, 10.77),  # TP.HCM
], crs='EPSG:4326')
print(f"Kiểu dữ liệu GeoSeries:\n{geo_series}\n")

Kiểu dữ liệu GeoSeries:
0    POINT (105.85 21.02)
1    POINT (106.63 10.77)
dtype: geometry

# Dữ liệu thành phố lớn Việt Nam (tọa độ mang tính tương đối). Dữ liệu này chỉ mang tính minh họa và nên được kiểm tra lại với dữ liệu chính thức.
vietnam_cities_data = {
    'city': ['Hà Nội', 'TP.HCM', 'Hải Phòng', 'Đà Nẵng', 'Cần Thơ', 'Biên Hòa', 'Huế', 'Nha Trang', 'Buôn Ma Thuột', 'Quy Nhon'],
    'province': ['Hà Nội', 'TP.HCM', 'Hải Phòng', 'Đà Nẵng', 'Cần Thơ', 'Đồng Nai', 'Thừa Thiên Huế', 'Khánh Hòa', 'Đắk Lắk', 'Bình Định'],
    'region': ['Miền Bắc', 'Miền Nam', 'Miền Bắc', 'Miền Trung', 'Miền Nam', 'Miền Nam', 'Miền Trung', 'Miền Trung', 'Miền Trung', 'Miền Trung'],
    'population': [8246600, 8993082, 2028514, 1134310, 1282937, 1104800, 455230, 423000, 340000, 284000],
    'longitude': [105.8542, 106.6297, 106.6881, 108.2022, 105.7469, 106.8439, 107.5905, 109.1967, 108.0373, 109.2189],
    'latitude': [21.0285, 10.8231, 20.8449, 16.0544, 10.0452, 10.9460, 16.4637, 12.2585, 12.6667, 13.7830],
    'is_port_city': [False, True, True, True, True, False, False, True, False, True],
}

# Tạo DataFrame thường trước
df = pd.DataFrame(vietnam_cities_data)
# tạo cột geometry từ longitude và latitude
geometry = [Point(xy) for xy in zip(df.longitude, df.latitude)]
# Tạo GeoDataFrame từ DataFrame và cột geometry
gdf = gpd.GeoDataFrame(
    df.drop(['longitude', 'latitude'], axis=1), 
    geometry=geometry, 
    crs='EPSG:4326'
)
gdf.head()

# Tạo GeoDataFrames từ polygons sử dụng shapely và gán vào geometry
# Tạo một số polygons ví dụ
geometry = [
    Polygon([(105, 20), (106, 20), (106, 21), (105, 21)]),
    Polygon([(106, 10), (107, 10), (107, 11), (106, 11)]),
    Polygon([(108, 16), (109, 16), (109, 17), (108, 17)]),
    Polygon([(109, 12), (110, 12), (110, 13), (109, 13)]),
    Polygon([(105, 10), (106, 10), (106, 11), (105, 11)]),
    Polygon([(106, 10), (107, 10), (107, 11), (106, 11)]),
    Polygon([(107, 16), (108, 16), (108, 17), (107, 17)]),
    Polygon([(109, 12), (110, 12), (110, 13), (109, 13)]),
    Polygon([(108, 12), (109, 12), (109, 13), (108, 13)]),
    Polygon([(109, 13), (110, 13), (110, 14), (109, 14)]),
]
polygon = gpd.GeoDataFrame(
    geometry=geometry,
    crs='EPSG:4326'
)
# Ta có thêm cột thuộc tính vào GeoDataFrame polygon
polygon["landcover"] = ["Urban", "Urban", "Rural", "Rural", "Urban", "Urban", "Rural", "Rural", "Rural", "Rural"]
polygon.head()

# Đọc dữ liệu từ local machine
districts = gpd.read_file(r'G:\My Drive\python\geocourse\data\vector\Vietnam_districts.geojson')
# Chuyển crs từ 4326 sang 32648 (UTM 48N)
districts = districts.to_crs(epsg=32648)
# Thêm cột diện tích 
districts['area_km2'] = districts.geometry.area/1e6  # Chuyển từ m2 sang km2
districts.head()

# Đọc dữ liệu từ url 
vietnam_data = gpd.read_file('https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_VNM_1.json')
vietnam_data.head()

# Lưu dữ liệu ra file GeoJSON
vietnam_data.to_file(os.path.join(outpath, 'vietnam_provinces.geojson'), driver='GeoJSON')

# Lưu dữ liệu ra file Shapefile
vietnam_data.to_file(os.path.join(outpath, 'Vietnam_provincess.shp'), driver='ESRI Shapefile')

# Lưu dữ liệu ra file Parquet
vietnam_data.to_file(os.path.join(outpath, 'Vietnam_provinces.parquet'), driver='Parquet')

# Lưu dữ liệu ra file GeoPackage
vietnam_data.to_file(os.path.join(outpath, 'Vietnam_provinces.gpkg'), driver='GPKG')

# Đọc dữ liệu từ url 
province = gpd.read_file('https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_VNM_1.json')
# Select two columns 
province = province[['VARNAME_1', 'geometry']]
# Rename the column
province = province.rename(columns={'VARNAME_1': 'province'})
province.head()

# Chuyển hệ tọa độ từ WGS84 (EPSG:4326) sang UTM 48N (EPSG:32648)
province_utm = province.to_crs(epsg=32648)
# Chọn tỉnh vĩnh phúc
vinhphuc = province_utm[province_utm['province'] == 'VinhPhuc']
# Tạo 1km buffer quanh tỉnh vĩnh phúc
vinhphuc_buffer = vinhphuc.buffer(1000)  # Buffer 1000 mét (1 km)
# Chuyển từ GeoSeries sang GeoDataFrame để dễ dàng xử lý
vinhphuc_buffer_gdf = gpd.GeoDataFrame(geometry=vinhphuc_buffer, crs='EPSG:32648')
vinhphuc_buffer_gdf.head()

# Đọc dữ liệu districts từ url
districts = gpd.read_file('https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_VNM_2.json')
# Đảm bảo dữ liệu vinhphuc có cùng hệ tọa độ với districts trước khi join
vinhphuc = vinhphuc.to_crs(districts.crs)
# Join dữ liệu tỉnh vĩnh phúc với dữ liệu districts để lấy ra các huyện thuộc tỉnh vĩnh phúc
vinhphuc_districts = gpd.sjoin(vinhphuc, districts, how='inner', predicate='intersects') # ngoài intersects còn có within, contains, touches, crosses, covers, covered_by.
vinhphuc_districts.head()