GeoPy1.1 地理数据处理入门

news2024/10/5 21:18:55

原作者：Damon 高校教师，中科院 GIS 博士

本文为原文章基础上，加上自己以及GPT4.0的总结整理而来

原活动链接

前言

在地理学及相关专业，我们往往需要处理大量的地球观测数据和模拟结果，基于这些数据的处理结果进一步分析，得到一些客观规律或者研究结论。但是在处理过程中，我们会发现这些数据可能来源不同，格式不一。对于一些常见的地理数据格式，比如栅格数据、矢量数据等，大家比较容易想到用专业软件来处理，比如 Arcgis、MapGIS、ENVI 等，这些软件使用简单，也很强大，但是也有着一定的不足，比如（1）往往只能处理专业格式的数据（2）批量处理一般有着一定的门槛，需要进行二次开发（3）一般为商业软件，正版的授权费用较高等等。

为了更方便地批量处理地理数据、与其他领域的数据进行联合分析、展示结果等，相信大家都会想到用编程来解决相关问题，而目前最火的编程语言就是 Python。Python 简单易学，尤其是对于地理数据处理来说，Python 中着大量的地学相关库，包括数据处理、可视化和模型开发等，方便我们在同一个工作环境下处理数据、分析数据、展示分析结果，研究结果也容易被别人复现，且方便我们在研究中使用最新的算法，比如深度学习等。如果你是一名正在或将要找工作的同学，那更好，学好 Python，offer 多多！

学习 Python 的资源在网上多如牛毛，然而这对于没有接触过编程的同学来说不一定是件好事儿，毕竟筛选资源也是一门技术活。其实，从实用性的角度来说，可以先学习能让你 “读懂” 别人代码的基础知识，然后剩下的高阶知识在不断的实践过程中去补充。

本关卡就是作为一个 “砖头”，学完本关，你可以做到从一堆杂七杂八的数据中筛选出 shp 格式的矢量数据，帮你完成数据的分类整理，踏出地理数据分析第一步。

import numpy as np
import pandas as pd

小练习：求一周的平均温度

假设某地一周的平均温度分别为 15.3、16.7、14.8、15.0、19.7、16.4、18.5 摄氏度，请计算本周的平均温度。

# a1至a7依次为周一到周日的平均温度
a1,a2,a3,a4,a5,a6,a7 = 15.3,16.7,14.8,15.0,19.7,16.4,18.5

mean_value = np.array([a1,a2,a3,a4,a5,a6,a7]).mean()  # 请在=后的空格处填入你的代码
print(mean_value)

16.62857142857143

通过赋值语句，可以将数据存储在变量中，并在程序中进行处理和操作，这是 Python 中常用的基本语法之一。

除了上面的简单数据，我们还可以把更为复杂的数据（比如一些地理数据）赋值给一个变量，比如下面的代码我们把一个栅格数据据赋值给了左边的变量 rs：

import rasterio as rio
from rasterio.plot import show
rs = rio.open('./data/GeoPy1/temp_grid_84/tem2015_84.tif')
show(rs);  # 仅展示rs变量中的数值数据

在这里插入图片描述

小练习：将文件夹下的文件路径都打印出来：

import glob

files = glob.glob('./data/GeoPy1/*')  # 获取文件夹下的所有文件
for i in files:
    print(i)

./data/GeoPy1\beijing.cpg
./data/GeoPy1\beijing.dbf
./data/GeoPy1\beijing.prj
./data/GeoPy1\beijing.qmd
./data/GeoPy1\beijing.shp
./data/GeoPy1\beijing.shx
./data/GeoPy1\beijing_temp.csv
./data/GeoPy1\hunan.cpg
./data/GeoPy1\hunan.dbf
./data/GeoPy1\hunan.prj
./data/GeoPy1\hunan.qmd
./data/GeoPy1\hunan.shp
./data/GeoPy1\hunan.shx
./data/GeoPy1\LC81970242014109LGN00.tif
./data/GeoPy1\Province_wgs84.cpg
./data/GeoPy1\Province_wgs84.dbf
./data/GeoPy1\Province_wgs84.prj
./data/GeoPy1\Province_wgs84.sbn
./data/GeoPy1\Province_wgs84.sbx
./data/GeoPy1\Province_wgs84.shp
./data/GeoPy1\Province_wgs84.shp.xml
./data/GeoPy1\Province_wgs84.shx
./data/GeoPy1\temp
./data/GeoPy1\temp_grid_84

🐋：glob 模块在后面关卡批量处理时会再次用到。

小练习：判断矢量和栅格并重命名

我们再来看另外一个案例：通过一个函数来给文件名字后面加上 “矢量” 或者 “栅格” 的标识，让业外人士也能清晰地分辨数据是矢量还是栅格。

def rename_file(x):
    """
    （这段话是块注释）
    根据x的名字来判断是矢量还是栅格
    如果是矢量，在文件名最前面加上"矢量"二字，
    如果是栅格，在文件名最前面加上"栅格"二字。
    """
    if x[-4:] == '.shp':
        x = '矢量' + x
    elif  x[-4:] == '.tif':
        x = '栅格' + x
    return x

# 调用函数并打印返回值
test_result = rename_file('china.shp')
print("新的文件名字为:", test_result)

新的文件名字为: 矢量china.shp

上面定义了一个很实用的函数，也许在工作和科研中会经常遇到，尤其是做数据入库的时候。

闯关题

按照题目要求完成以下任务

写一个函数，输入 “china.shp”，然后函数进行判断该文件是否是矢量数据，如果是矢量数据（假设是 ESRI shapefile），返回该文件名字（包括后缀），如果不是，则返回 “文件输入错误”，然后将函数的返回值赋值给变量 a1。

可以通过判断文件后缀来确定是否是矢量文件（比如是否等于 “.shp”）
split() 函数可以对 str（字符串）对象进行分割
split() 函数分割字符串之后得到一个 list，可以通过下标（也叫索引）来 “取出” 对应位置的元素
list 的下标是从 0 开始

# ...你的代码...
def is_shp(x):
    file_type = x.split('.')[1] # 请将'_'替换为合适的下标数字
    if file_type == 'shp':
        return x
    else:
        return '文件输入错误'

a1 = is_shp('china.shp')  # 请将'_'替换为给定的文件名
a1