मेरे पास वैश्विक समुद्र विज्ञान (ओमेगा) डेटा की एक नेटसीडीएफ फ़ाइल है जो अपेक्षाकृत मोटे स्थानिक रिज़ॉल्यूशन पर 33 गहराई स्तरों के साथ है। मेरे पास बहुत बेहतर रिज़ॉल्यूशन पर एक वैश्विक बाथमीट्री रेखापुंज भी है। मेरा लक्ष्य नेटसीडीएफ फ़ाइल से सीबेड ओमेगा डेटा प्राप्त करना है, वांछित गहराई निर्धारित करने के लिए बाथमीट्री डेटा का उपयोग करना है। मेरा कोड अब तक;

    library(raster)
    library(rgdal)
    library(ncdf4)

# Aragonite data. Defaults to CRS WGS84
ncin <- nc_open("C:/..../GLODAPv2.2016b.OmegaA.nc")
ncin.depth <- ncvar_get(ncin, "Depth")# 33 depth levels  

omegaA.brk  <- brick("C:/.../GLODAPv2.2016b.OmegaA.nc")
omegaA.brk  <-rotate(omegaA.bkr)# because netCDF is in Lon 0-360. 

# depth raster.  CRS WGS84
r<-raster("C:/....GEBCO.tif") 

# resample the raster brick to the resolution that matches the bathymetry raster
omegaA.brk  <-resample(omegaA.brk, r, method="bilinear")

# create blank final raster 
omegaA.rast         <- raster(ncol = r@ncols, nrow = r@nrows)
extent(omegaA.rast) <- extent(r)
omegaA.rast[]       <- NA_real_

#  create vector of indices of desired depth values
depth.values<-getValues(r)
depth.values.index<-which(!is.na(depth.values))


# loop to find appropriate raster brick layer, and extract the value at the desired index, and insert into blank raster

for (p in depth.values.index) { 
  dep.index         <-which(abs(ncin.depth+depth.values[p]) == min(abs(ncin.depth+depth.values[p]))) ## this sometimes results in multiple levels being selected

  brk.level         <-omegaA.brk[[dep.index]]  # can be more than on level if multiple layers selected above. 
  omegaA.rast[p]    <-omegaA.brk[[1]][p] ## here I choose the first level if multiple levels have been selected above

  print(paste(p, "of", length(depth.values.index))) # counter to look at progress. 
}

समस्या: परिणाम बड़े पैमाने पर अंतराल (एनए) के साथ एक रेखापुंज है जहां डेटा होना चाहिए। अंतराल अक्सर एक विशिष्ट आकार लेते हैं - उदाहरण के लिए, एक समोच्च का पालन करें, या एक लंबी सीधी रेखा के साथ। मैंने एक फसली उदाहरण चिपकाया है।

छवि विवरण यहां दर्ज करें

मुझे लगता है कि ऐसा इसलिए हो सकता है क्योंकि या तो 1) किसी कारण से लूप में 'कौन सा' कथन मेल नहीं ढूंढ रहा है या 2) अनुमानों का एक गलत संरेखण बनाया गया है जिसे मैंने पढ़ा है 'रोटेट' का उपयोग करते समय हो सकता है।

मैंने यह सुनिश्चित करने की कोशिश की है कि सभी विस्तार, संकल्प, कोशिकाओं की संख्या, और सीआरएस सभी समान हैं, जो वे प्रतीत होते हैं।

प्रक्रिया को तेज करने के लिए मैंने अपनी रुचि के क्षेत्र में वैश्विक ईंट और बाथ रैस्टर को क्रॉप किया है, फिर से जाँच कर रहा हूँ कि सभी स्थानिक संकल्प, आदि आदि मेल खाते हैं - मैंने उन चरणों को यहाँ सरलता के लिए शामिल नहीं किया है।

नुकसान में। किसी भी मदद का स्वागत है!

1
user2175481 4 जिंदा 2018, 14:52

1 उत्तर

सबसे बढ़िया उत्तर

एक प्रतिलिपि प्रस्तुत करने योग्य उदाहरण के बिना, इस तरह की समस्याओं को हल करना मुश्किल है। मैं यह नहीं बता सकता कि आपकी समस्या कहाँ है, लेकिन मैं आपके सामने वह दृष्टिकोण प्रस्तुत करूँगा जो मैं कोशिश करूँगा। शायद यह अच्छा है, शायद यह बुरा है, मुझे नहीं पता लेकिन यह आपको अपनी समस्या के समाधान का रास्ता खोजने के लिए प्रेरित कर सकता है।

मेरी समझ से, आपके पास ओमेगा (33 परतें/गहराई) की एक ईंट और एक बाथमीट्री रेखापुंज है। आप समुद्र के तल पर ओमेगा मान प्राप्त करना चाहते हैं। यहां बताया गया है कि मैं कैसे करूंगा:

  1. ओमेगा को एक ही रिज़ॉल्यूशन और बाथमेट्री के लिए एक रेखापुंज बनाएं
  2. बाथिमेट्री रेखापुंज को 0-1 की 33 तीन परतों वाली रेखापुंज ईंट में रूपांतरित करें। जैसे यदि एक विशेष पिक्सेल के लिए समुद्र तल 200 मीटर पर है, तो 200 के अलावा अन्य सभी गहराई परत पर यह पिक्सेल 0 है और 200 के लिए 1 है। इसे प्रोग्राम करने के लिए, मैं लंबा रास्ता तय करूंगा, कुछ ऐसा

:

r_1 <- r
values(r_1) <- values(r)==10  # where 10 is the depth (it could be a range with < or >)
r_2 <- r
values(r_2) <- values(r)==20
...
r_33 <- r
values(r_33) <- values(r)==250
r_brick <- brick(r_1, r_2, ..., r_33)
  1. फिर आप अपने दोनों रेखापुंज ईंटों को कई गुना करते हैं। उनका एक ही आयाम है, यह आसान होना चाहिए। आउटपुट 33 परतों की एक रेखापुंज ईंट होनी चाहिए जिसमें 0 हर जगह 0 के साथ हो, जहां यह समुद्र के नीचे नहीं है और कहीं और ओमेगा का मूल्य नहीं है।
  2. पहले प्राप्त ईंट की सभी परतों को एक योग के साथ एक साधारण रेखापुंज में मिलाएं।

यह काम करना चाहिए। यदि आपको रेखापुंज ईंट से निपटने में समस्या है, तो आप डेटा को आधार R सरणियों में बना सकते हैं, यह सरल हो सकता है।

आपको कामयाबी मिले।

0
Bastien 4 जिंदा 2018, 17:31