xuefliang

import pandas as pd data = pd.DataFrame({     "Student": ["A", "B", "C", "D"],      "Physics": [89,34,23,56],      "Chemistry": [34,56,98,56],      "Math": [34,94,50,59],     'Name':["LI","Zhang","Wang","Xiao"]     }) pd.eval("data.Physics+data.Chemistry+data.Math") data.eval("total=Physics+Chemistry+Math").query('Math>50') data.eval('total2=Student',engine='python') data.eval('total2=Student+Name',engine='python') # 'numexpr': This default engine evaluates pandas objects using # numexpr for large speed ups in complex expressions with large frames. # 'python': Performs operations as if you had eval’d in top # level python. This engine is generally not that useful. (data.eval('total=sqrt(Math)')  .query('Chemistry>34')) data=pd.read_excel(r"C:\Users\xuefe\Desk

 from sklearn.preprocessing import LabelEncoder #Auto encodes any dataframe column of type category or object. def dummyEncode(df):         columnsToEncode = list(df.select_dtypes(include=['category','object']))         le = LabelEncoder()         for feature in columnsToEncode:             try:                 df[feature] = le.fit_transform(df[feature])             except:                 print('Error encoding '+feature)         return df      X=dummyEncode(X)

 import cx_Oracle import os import pandas as pd import numpy as np import geopandas as gpd yg=pd.read_excel(r"C:\Users\xuefe\Documents\乙肝项目.xlsx") yg.info() # duplicated 判断是否有重复项 mask=yg[['GRDA_CODE','YM_MC']].duplicated(keep=False) #选择重复 yg[mask] yg.duplicated() #去重 yg[~mask] yg['YM_MC'].drop_duplicates() yg.drop_duplicates(['YM_MC']) # 当keep=False时，就是去掉所有的重复行  # 当keep=‘first'时，就是保留第一次出现的重复行  # 当keep='last'时就是保留最后一次出现的重复行。

 sqltxt = '''   select t.grda_code,        t.grda_xm,        e.csrq,        t.ym_mc,        t.jz_sj,        t.jz_zc,        t.ym_yxq,        t.ym_ph,        t.ym_bm,        t.ymsccj_mc,        t.jzdd_mc,        t.jzdd_dm   from inoc_jzjl t, grda_cr e  where e.grda_cr_lsh = t.grda_et_lsh  and t.config_ymxx_lsh in ('2901','2902','2903','2904','2905','2906')  and  t.jz_sj>=to_date('2021-01-01','yyyy-mm-dd')  and jcqk='1'  '''   jzjl=pd.read_sql(sqltxt,con) jzjl.columns=[str(item).lower() for item in jzjl.columns] jzjl['jz_sj']=jzjl['jz_sj'].dt.date jzjl['csrq']=pd.to_datetime(jzjl['csrq'],format='%Y%m%d') jzjl=jzjl[jzjl['ym_ph'].isin(['20210411','20210310','20210205','20210206','20201212','20210102'])] test=jzjl test['pc']=np.where(test['ym_ph'].isin(['20210411','202103

  conda install contextily --channel conda-forge

 import cx_Oracle import os import re import math import pandas as pd import numpy as np from sqlalchemy import create_engine from sqlalchemy.engine import url from matplotlib import pyplot as plt  os.environ['NLS_LANG'] = 'AMERICAN_AMERICA.ZHS16GBK' con = cx_Oracle.connect("ipvsdb", "Mygh@2021#", "192.168.30.48/JZDB1") sqltxt = '''  select i.grda_code,        e.zjhm,        e.csrq,        i.jz_zc,        i.jz_sj,        i.jzdd_dm,        i.ym_bm   from inoc_jzjl i, grda_cr e  where e.grda_cr_lsh = i.grda_et_lsh    and e.jcqk = '1'    and i.config_ymxx_lsh in ('1801','1802','1803')    and i.jz_sj >=to_date('2020-01-01','yyyy-mm-dd')  '''   cxr=pd.read_sql(sqltxt,con) cxr.columns=[str.lower(item) for item in cxr.columns] cxr['shi']=cxr['jzdd_dm'].str[0:4] cxr['xian']=cxr['jzdd_dm'].str[0:6] cxr['age']=((cxr['jz_sj']-cxr[

 library(sf) library(tidyverse) library(magrittr) `%nin%` = Negate(`%in%`) xz <- st_read("D:\\乡镇边界\\乡镇边界.shp") %>%    select(Name) %>%    st_make_valid() %>%  #几何体修正几何图形   st_write("D:\\乡镇边界\\乡镇边界.geojson") xz %>%    st_make_valid()  xz %>%    st_is_valid() %>%    as_tibble() -> tst xz %>%    filter(na.omit(st_is_valid(xz)) == TRUE) -> xz #过滤冒犯的几何图形(通过基于结果的空间对象的子集) # empty geometries, using any(is.na(st_dimension(x))) # corrupt geometries, using any(is.na(st_is_valid(x))) # invalid geometries, using any(na.omit(st_is_valid(x)) == FALSE); in case of corrupt and/or invalid geometries, # in case of invalid geometries, query the reason for invalidity by st_is_valid(x, reason = TRUE) # you may be succesful in making geometries valid using st_make_valid(x) or, if st_make_valid is not supported by # st_buffer(x, 0.0) on non-corrupt geometries (but beware of the bowtie example above, where st_buffer removes one half). # After succesful a st_mak

 library(sf) library(tidyverse) library(magrittr) `%nin%` = Negate(`%in%`) ts <- st_read("C:\\Users\\ruiying\\Desktop\\甘肃乡镇\\gansu.TAB") %>%    filter(str_detect(地区编码,"^620502") | str_detect(地区编码,"^620503")) %>%   rename(code=地区编码) %>%    filter(row_number() %nin% c(44,45)) # ts %>%  #   st_transform(4326) # ts <- st_read("C:\\Users\\ruiying\\Desktop\\ts.geojson") # bl <- read.csv("C:\\Users\\ruiying\\Desktop\\报告卡2021-11-03（8时）.csv") # bl %<>% #   mutate(xiang=str_sub(bl$现住详细地址,10,12)) %>% #   group_by(xiang) %>% #   summarise(n=n()) xiang <- read.csv("C:\\Users\\ruiying\\Desktop\\xiang.csv") %>%   select(-X) ts %>%    left_join(xiang,by=c("NAME"="NAME")) %>%    ggplot() +   geom_sf(aes(fill=n))+   scale_fill_gradientn(breaks = c(2,4,6),                        colors=c("gray", "yellow", "orange", "red"))+   coord_sf(xlim =

 library(tidyverse) dat = data.frame(y1 = rnorm(10),y2 = rnorm(10)+10) re = list() #1 for(i in 1:2){ #2   re[[i]] = mean(dat[,i]) #3 }  map(dat,mean) map_dbl(dat,mean) map_df(dat,mean) map_df(dat,~mean(.)) iris %>%    split(.$Species) %>%    map(.,~aov(Sepal.Length ~ Sepal.Width,data=.) %>% summary) # aov(Sepal.Length ~ Sepal.Width,data = iris) mtcars %>%    split(.$cyl) %>%    map(.,~lm(mpg~wt+hp+drat,data = .)) %>%    map(coef) %>%    map_df('wt') mtcars %>%    split(.$cyl) %>%    map(.,~lm(mpg~wt+hp+drat,data = .)) %>%    map(coef) %>%    as.data.frame() x <- list(   list(-1, x = 1, y = c(2), z = "a"),   list(-2, x = 4, y = c(5, 6), z = "b"),   list(-3, x = 8, y = c(9, 10, 11)) ) # select by name map_dbl(x,'x') # select by position map_dbl(x,1) # select by name and position map_dbl(x,list('y',1)) wt <- c(5,  5,  4,  1)/15 x <- c(3.7,3.3,3.5,2.8) weighted.mean(x, wt)  # 3.7*0.3333+3.3*0.3333+3.5*0.2

 library(sf) library(tidyverse) 读写geojson shi <- st_read("C:\\Users\\ruiying\\Desktop\\shi.geojson") st_write(shi, "C:\\Users\\ruiying\\Desktop\\shi2.geojson",append=F) 读写json weihai <- st_read("D:\\使用GEOJSON数据绘制地图\\威海市.json") st_write(weihai %>% select(name,adcode),"C:\\Users\\ruiying\\Desktop\\weihai.geojson") st_write(meuse_sf,"weihai.shp") st_write(meuse_sf,"weihai.geojson",append = F) st_write(meuse_sf,"weihai.json",driver = 'GeoJSON',append = F) st_write(meuse_sf,"weihai.gpkg",append = F) test <- st_read('C:\\Users\\ruiying\\Desktop\\shi.geojson')

 shi %>%    left_join(shi_df,by=c('市代码'='市代码')) %>%    filter(str_detect(市代码,'^62')) %>%    ggplot()+   geom_sf(aes(fill=bl),size=0.01)+   scale_fill_gradientn(breaks = c(15,30,45),                        colors=c("gray", "yellow", "orange", "red"))

 library(sampling) mtcars_order <- mtcars %>% arrange(gear) mtcars_order %>% group_by(gear) %>% summarise(n=n()) smn <- round(prop.table(table(mtcars_order$gear))*10) stratID <- strata(mtcars_order,stratanames='gear',size = smn,method = 'srswor') head(stratID) strat<-getdata(mtcars_order,stratID) library(doBy) mtcars %>% group_by(am,gear) %>% summarise(n=n()) mtcars %>% group_by(am) %>% summarise(n=n()) sampleBy(formula = ~ iris$Species,frac=0.1,data=iris) sample_by(mtcars,formula = ~am+gear,frac = 0.2)

 df <- data.frame(id = 1:15,                  grp = rep(1:3,each = 5),                   frq = rep(c(3,2,4), each = 5)) set.seed(22) df %>%   group_by(grp) %>%                                 # for every group   summarise(d = list(data.frame(id=id)),            # create a data frame of ids             frq = unique(frq)) %>%                  # get the unique frq value   mutate(v = map2(d, frq, ~sample_n(.x, .y))) %>%   # sample using data frame of ids and frq value   unnest(v) %>%                                     # unnest sampled values   select(-frq)                                      # remove frq column (if needed)

 xzqh <- read.csv('./Documents/地区编码.csv', fileEncoding = "GBK",stringsAsFactors=FALSE) # locale:locale参数是 readr 包中很重要的一个参数，指定日期使用的月和日的名称，时区，字符编码，日期格式，数字的小数和点位数和分隔符。 xzqh <- read_csv('./Documents/地区编码.csv',                   locale = locale(date_format = "%Y/%m/%d",                                  time_format = "h%Hm%Ms%S %p",                                  encoding = 'GBK',                                  tz = 'Asia/Shanghai'))

  import  numpy  as  np import  pandas  as  pd df=pd.read_excel( r 'C: \U sers \r uiying \D esktop \评 审结果 . xlsx' ) df[ 'xiang' ]=df[ 'jgbm' ].str[ 0 : 8 ] df[ 'xiang' ].count() df[ 'xiang' ].value_counts() df[df[ 'pdjb_bm' ]>= 3 ][ 'pdjb_bm' ].count() df[ 'xiang' ].drop_duplicates().count() df.describe() result=df.groupby( 'xiang' ).apply( lambda  x: x[x[ 'pdjb_bm' ]>= 3 ][ 'pdjb_bm' ].count()).to_frame() type (result) result.head() result.columns=[ 'xiang' , 'num' ] result.describe() result.columns result.rename( columns ={ 'xiang' : 'xian' }, inplace = True ) result.to_csv( r 'C: \U sers \r uiying \D esktop \结 果 . csv' )

  import  cx_Oracle import  os import  re import  math import  pandas  as  pd import  numpy  as  np from  sqlalchemy  import  create_engine from  sqlalchemy.engine  import  url from  id_validator  import  validator from  matplotlib  import  pyplot  as  plt  os.environ[ 'NLS_LANG' ] =  'AMERICAN_AMERICA.ZHS16GBK' sqltxt= ''' select jz_sj,ym_mc,jzdd_dm from inoc_jzjl where jz_sj >=to_date('2020-01-01','yyyy-mm-dd') and jz_sj <=to_date('2020-08-30','yyyy-mm-dd') and sfmf=1 ''' con = cx_Oracle.connect( "****" ,  "********" ,  "192.168.30.48/JZDB1" ) jzjl = pd.read_sql(sqltxt,con)

import pandas as pd import numpy as np from sklearn import datasets, preprocessing from sklearn.model_selection import train_test_split #%% iris = datasets.load_iris() x = pd.DataFrame(iris['data'], columns=iris['feature_names']) y = pd.DataFrame(iris['target'].astype(str), columns=['target_names']) df = pd.concat([x,y], axis=1) #%% '檢查&補值' df.isna().sum() df=df.fillna(0) #0補值 df=df.fillna(method='pad') #前一筆 df=df.fillna(method='bfill') #後一筆 df['target_names']=df['target_names'].fillna((df['target_names'].mode())) #眾數 df['sepal length (cm)']=df['sepal length (cm)'].fillna((df['sepal length (cm)'].mean())) #平均值 #%% '資料轉換-標準化' min_max_scaler = preprocessing.MinMaxScaler() df = min_max_scaler.fit_transform(df) #%% '單熱' one_hot = pd.get_dummies(df) labels = np.array(one_hot['sepal length (cm)']) #應變數 features= one_hot.drop(&

 import  pandas as pd import numpy as np df=pd.DataFrame({'amount':[100,200,300,400,500],'list':['','商品1','商品2','','商品3']}) df['x1'] = df.apply(lambda x: x.amount if x.list != "" else 0, axis=1) df['x2']=np.where(df['list']=='',0,df['amount']) df.loc[df['list']!='','x3']=df['amount'] df=df.fillna(0)

 cr %>%    filter(str_detect(shi,'^62')) %>%    group_by(shi) %>%    summarise(v1=sum(year>=18,na.rm = T),v2=sum(year<18,na.rm = T),v3=sum(year>=0),n=n()) cr %>%    mutate(age_group=case_when(     year<6 ~ '<6',     year>=6 & year<12 ~ '6-11',     year>=12 & year<15 ~ '12~14',     year>=15 & year<18 ~ '15-17',     year>=18 & year<59 ~ '18-59',     year>=60 & year<70 ~ '60-69',     year>=70 ~ '>=70')) %>%    filter(str_detect(shi,'^62')) %>%    group_by(age_group) %>%   summarise(n=n())   

 first_category_name = list.files("./Documents/礼县附件二/")             #list.files命令得到文件夹下所有文件夹的名称【修改】 dir = paste("./Documents/礼县附件二/",first_category_name,sep="")    #用paste命令构建路径变量dir, n = length(dir)    datalist = list() for(i in 1:n){   if(str_detect(dir[i],'xlsx')){     datalist[[i]]=readxl::read_xlsx(dir[i],skip = 3)     }   else{     datalist[[i]]=readxl::read_xls(dir[i],skip = 3)   } } nrows=sapply(datalist, nrow) sum(nrows) n = length(dir)    i=17 datalist = data.frame() if (str_detect(dir[i], 'xlsx')) {   datalist = readxl::read_xlsx(dir[i], skip = 2, sheet = 1) } else {   datalist = readxl::read_xls(dir[i], skip = 2, sheet = 1) }

 data(iris) iris <- iris[,-c(5)] library(rminer)  model <- fit(Sepal.Length~.,data=iris,model="svm",              kpar=list(sigma=0.1),C=1) svm <- Importance(model,iris,measure="AAD") svm$imp par(mfrow=c(1,3)) boxplot(iris$Sepal.Width)$out boxplot(iris$Petal.Length)$out boxplot(iris$Petal.Width)$out #因為離群值不多，所以先不轉 n <- nrow(iris) set.seed(1117) subiris <- sample(seq_len(n), size = round(0.7 * n)) traindata <- iris[subiris,] testdata <- iris[ - subiris,] library(e1071) tune.model <- tune.svm(Sepal.Length~.,data=traindata,type="eps-regression",kernel="radial",                        range=list(cost = 2^c(-8,-4,-2,0), epsilon = seq(0,10,0.1))) tune.model$best.model #挑 lowest MSE  model <- svm(Sepal.Length~.,data=traindata,type="eps-regression",kernel="radial",cost=1,epsilon=0.1,gamma=0.3333333) future <- predict(model,testdata) future <- as.data.frame(future) final <- cbind(future,testdata) lib

 #GB library(gbm) set.seed(123) model <- gbm(formula = Sepal.Length ~ .,distribution = "gaussian",              data = iris,n.trees = 479,              interaction.depth = 1,              shrinkage = 0.1,n.minobsinnode = 15,              bag.fraction = 1,              train.fraction = 0.5,cv.folds = 10) #分類用 bernoulli 迴歸用gaussian(RMSE) #XGB data("iris") par(mfrow=c(1,4)) boxplot(iris$Sepal.Length)$out boxplot(iris$Sepal.Width)$out boxplot(iris$Petal.Length)$out boxplot(iris$Petal.Width)$out #因為離群值不多，所以先不轉 n <- nrow(iris) set.seed(1117) subiris <- sample(seq_len(n), size = round(0.7 * n)) traindata <- iris[subiris,] testdata <- iris[ - subiris,] features <- setdiff(names(iris), "Sepal.Length") library(vtreat)#單熱  library(dplyr) onehot <- vtreat::designTreatmentsZ(iris, features, verbose = FALSE) 單熱 <- onehot %>%   magrittr::use_series(scoreFrame) %>%           dplyr::filter(code %in% c("clean", "lev")) %>

 library(naniar)  data("iris") iris <- iris[,-c(5)] any_na(iris) library(randomForest) model <- randomForest(Sepal.Length ~ Sepal.Width+Petal.Length+Petal.Width, data = iris,                          importane = T, proximity = T,                        do.trace = 100,na.action = na.roughfix) par(mfrow=c(1,1)) plot(model)  importance(model) library(rminer) rf <- Importance(model,iris,measure="AAD") rf$imp par(mfrow=c(1,3)) boxplot(iris$Sepal.Width)$out boxplot(iris$Petal.Length)$out boxplot(iris$Petal.Width)$out #因為離群值不多，所以先不轉 n <- nrow(iris) set.seed(1117) subiris <- sample(seq_len(n), size = round(0.7 * n)) traindata <- iris[subiris,] testdata <- iris[ - subiris,] features <- setdiff(x = names(traindata), y = "Sepal.Length") par(mfrow=c(1,1)) set.seed(123) tuneRF(x = traindata[features], y = traindata$Sepal.Length,        mtryStart = 1,ntreeTry = 450) set.seed(123) model <- randomForest(Sepal.Length ~ Sepal.Width+Petal.Length+Pe

library(naniar)  data("iris") iris <- subset(iris,Species!="setosa") any_na(iris) n <- nrow(iris) set.seed(1117) new <- iris[sample(n),] t_idx <- sample(seq_len(n), size = round(0.7 * n)) traindata <- iris[t_idx,] testdata <- iris[ - t_idx,] library(rpart)  library(rpart.plot) dtreeM <- rpart(Species ~ ., data = traindata,                  method = "class") #數值型改anova rpart.plot(dtreeM,digits=2,varlen=20)#圖 printcp(dtreeM) #最佳cp=min(xerror) 每次分割能改善模型 剪枝 <- prune(dtreeM,cp = dtreeM$cptable[which.min(dtreeM$cptable[,"xerror"]),"CP"]) future <- predict(剪枝 , testdata, type = "class") future <- as.data.frame(future) final <- cbind(future,testdata) confusion <- table(final$Species,final$future, dnn = c("实际", "预测")) confusion accuracy <- sum(diag(confusion)) / sum(confusion) accuracy  

 data(iris) iris <- iris[,-c(5)] #找自變數 library(rminer) library(tidyverse) model <- lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris) lm <- Importance(model,iris,measure="AAD") lm$imp #檢查離群值 par(mfrow=c(1,3)) boxplot(iris$Sepal.Width)$out boxplot(iris$Petal.Length)$out boxplot(iris$Petal.Width)$out par(mfrow=c(1,3)) scatter.smooth(x=iris$Sepal.Width, y=iris$Sepal.Length) scatter.smooth(x=iris$Petal.Length, y=iris$Sepal.Length) scatter.smooth(x=iris$Petal.Width, y=iris$Sepal.Length) #隨機抽樣 n <- nrow(iris) set.seed(1117) subiris <- sample(seq_len(n), size = round(0.7 * n)) traindata <- iris[subiris,] testdata <- iris[ - subiris,] #建模 残差的三个检验>0.05 model <- lm(Sepal.Length~poly(Sepal.Width,3)+Petal.Length+Petal.Width,traindata) #條件2~4 library(car) ncvTest(model)#>a 殘差變異數有同質性 shapiro.test(model$residuals) #>a 殘差常態 library(lmtest) dwtest(model)#>a 殘差獨立 vif(model)  #<10 ok 10~100可能過度配適 #預測 future <- predict(model,testdata) future <

 library(glmnet) library(tidyverse) data("iris") iris <- iris[,-c(5)] #自变量重要性 #ridge 自变量重要性变不为0 ,lasso 自变量重要性变为0 par(mfrow=c(1,3)) boxplot((iris$Sepal.Width))$out boxplot(iris$Petal.Length)$out boxplot(iris$Petal.Width)$out set.seed(1117) subiris <- sample(seq_len(nrow(iris)),size = round(0.7*nrow(iris))) traindata <- iris[subiris,] %>% as.matrix() testdata <- iris[-subiris,] %>% as.matrix() trainx <- traindata[,c(2:4)] trainy <- traindata[,c(1)] testx <- testdata[,c(2:4)] testy <- testdata[,c(1)] ridge <- cv.glmnet(x=trainx,y=trainy,alpha=0) #alpha=0为ridge,=1为lasso,k=10,交叉验证 #視覺化&選自變量 coef(ridge, s = "lambda.min") %>%   as.matrix() %>%   as.data.frame() %>%   add_rownames(var = "var") %>%   rename(coef=`1`) %>%    filter(var != "(Intercept)") %>% #剔除截距項   top_n(3, wt = coef) %>%   ggplot(aes(coef, reorder(var, coef))) +   geom_bar(stat = "identity", width=0.2,            color=&

  df03 <- billboard[,3] separate(data = df03,          col = date.entered,          into = c("year", "month", "day"),          sep = "\\-",          remove = F) -> df04 df04 unite(data = df04,       col = date,       year, month, day,       sep = "/") df05 <- read_xlsx('./Documents/批号价格.xlsx',sheet = 2) df05 %>%    fill(c(企业名称,批号,`价格/剂（元）`),.direction='down') %>%    unite(col = id,企业名称,批号,sep='-')

 HKEY_CURRENT_USER\SOFTWARE\Microsoft\Windows\CurrentVersion\Internet Settings Find a file with the name 'ProxyServer' and delete it.

 library(sf) library(gstat) library(openxlsx) library(magrittr) library(tidyverse) mycrs <- '+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +datum=WGS84 +no_defs' data_jiangshui <- read.xlsx('./Documents/wang_map/降水2005-2009修订后.xlsx') gs_station<- st_read("./Documents/wang_map/甘肃站点.shp")  %>%    st_transform(mycrs) gs_station %<>%    inner_join(data_jiangshui,by=c('stationid'='stationid'))  gs<- st_read("./Documents/wang_map/town_62.shp") %>%    st_transform(mycrs) gs.border <- st_boundary(gs) gs.bou <- st_union(gs) %>%    st_boundary() plot(gs.bou) grid <- st_make_grid(gs.bou, n = c(50, 50)) %>%    st_transform(st_crs(gs_station)) plot(grid) st_crs(gs_station)$proj4string st_crs(grid)$proj4string vg.gs <- variogram(gs_station$`2005-7` ~ 1, gs_station) # vg.gs # plot(vg.gs$dist, vg.gs$gamma) # vgm.gs <- vgm(80, "Gau", 500,nugget=200) fit.gs &

 library(automap) library(tidyverse) library(openxlsx) library(magrittr) library(raster) library(sf) mycrs <- '+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +datum=WGS84 +no_defs' gansu <- st_read("~/wang_map/甘肃站点.shp") jiangshui <- read.xlsx('~/wang_map/降水2005-2009修订后.xlsx') gansu %<>%    inner_join(jiangshui,by=c('stationid'='stationid')) %>%   st_transform(crs=mycrs) %>%    as_Spatial()  gsmap <- st_read("~/wang_map/town_62.shp") %>%    st_transform(crs=mycrs) %>%    as_Spatial()  #甘肃轮廓融合 # gsmapline <- rgeos::gUnaryUnion(gsmap) grd <- makegrid(gsmap, n = 10000) names(grd) <- c("X", "Y") gridded(grd)  <- ~ X+Y # Create SpatialPixel object plot(grd) # coordinates(grd) <- c("X", "Y") # gridded(grd) <- TRUE  # Create SpatialPixel object # fullgrid(grd)  ######### 泛克里金法grd需要指定 grd <- as.data.frame(spsample(