texts: Analyze relationship between two layers of points?

texts:
I have a layer of points that are patients who were readmitted to the hospital and a layer of points of patients who were not readmitted to the hospital.

I want to find geographical factors that influence this readmission. Some of the factors that I have GIS information: distance to fast food, fitness centers, churches, grocery stores

Using distance, what would be an efficient and logical way to analyze this data on QGIS?


texts:
I'm doing a lab and I'm trying to understand exactly why I'm using the tool "minus".

I have a raster which I'm calling "heights". I also have two rasters with water surfaces and watercourse. I have changed the values in the water rasters so that water pixels has the value 10 and noData has the value 0.

I have then used the tool "plus" to put these two rasters together into a joint raster which I'm calling "water" plain and simple.

Now, I am supposed to use the tool "Minus" to put in the water raster together with the raster heights, I can't understand why I need to do that.

My first thought was that the intention is to make sure that it's only the watersurfaces that shows up in my result later on and not land surfaces, but that doesn't seem right.

Why am I using minus in this scenario? It's to lower the watercourses I guess but why is that needed?

Update: The final goal of the lab is to find out from what area of the map a poison is coming from. To investigate this I have made an analysis of the watersurfaces in the area based on the water surfaces data and the flow direction.

Update: When answering the question "why do you have to use minus in this scenario " What I answered was basically

"after I have run the layer heights minus the layer water the watercourses will be lower than they really are in comparison to the surrounding surface. This happens because the pixels in the raster water i have the water values at 10 everything else at 0 (with reclassify). The height data will then get a value that's -10 lower than all the water courses. I'm doing this to get elevation in all the watercourses and in the modell the water will prefer to travel trough the water courses that's within the flow direction area. If I don't do this then some water courses may be without elevation and /or streams will be in the wrong places, even if I have a perfect height model, compared to where the actual streams are. So, I'm using this method to increase the chance that the watercourses are where they should be"

This was dead wrong. The answer I received that it will very well work with a perfect height model and that the height model should do just fine for analyzing. So the answer was wrong, that's not why I was using minus.

I'm totally lost here. No idea why of would use minus if that's not the answer. Any ideas?


texts:
I already have an ASP.net application using VS 2010 and .NET framework version 4 and ArcGis Server 10, and i need to add just a page, have map and toolbar.

When I try to add map control and run the page the map is not displayed. So i just want to add page to that application having map and toolbar controls.So can anyone can suggest a solution?

Thank you in advance.


texts: My Coordinate Reference System listing is in the wrong language

texts:
I just installed QGIS Desktop 2.0.1 on a desktop (Windows 8-64).

All the menus are in English-US (as they should be), but the CRS listings -- when I try to change them for a project are a layer -- are in French.  I fiddled with changing the language parameters in Settings | Options | Locale, but no matter what I do, the CRS listings still come up in French.

You can see in the attached image that the menus are in English, but the CRS names are in French (i.e., "Projection conique conforme" rather than "Conformal conic projection").

Any idea what I need to do to fix this?

Dave


texts: Reading GeoPackage file with GeoTools: Error on missing gpkg_contents database

texts:
I am trying to read a GeoPackage file using GeoTools with:

GeoPackage gp = new GeoPackage(new File("myFile.gpkg"));
FeatureEntry fe = gp.features().get(0);
SimpleFeatureReader fr = gp.reader(fe, filter, new DefaultTransaction());

while(fr.hasNext()) {
    SimpleFeature sf = fr.next();
    // Do things...
}
fr.close();
gp.close();

I get the following error:

  [SQLITE_ERROR] SQL error or missing database (no such table:
  gpkg_contents)

which I have no clue how to handle. The input file was produced from QGIS.
Is there something to take care of concerning this gpkg_contents database?


texts: Free heightmap dataset for Kraków, Poland?

texts:
Is there a freely available (governmental?) heightmap dataset for Kraków city (in Poland) with horizontal grid ~5m (or better) and vertical resolution ~ +/- 0.1m? Ideally, it would include building heights (i.e. not "ground level").

If you don't know one, at least do you know of some good open forum where I could ask such a question with some hope of getting an answer?

Sorry if that question doesn't use proper terms and/or is too unclear, or otherwise totally ridiculous; I'm totally unexperienced in GIS.


texts: Intersect tool in ArcGIS merges output points if they intersect same line segment

texts:
I wrote a quick script to automate the intersect geoprocessing. There are 2 inputs, Polyline A & Polyline B, and 1 output- a point file with XY where the two lines intersect.
Problem is if the fist input line (A) crosses the second input line(B) twice,
 it creates 2 points but "merges" them in the output table. The XY calculated thus, is halfway between these two points. What I am looking for is 2 seperate points in the output table with correct XY.
(I tried it with manual intersect tool from toolbox and it does the same thing)


texts: Rebuilding indexes in file geodatabase using ArcPy

texts:
I am trying to automate the rebuild of file geodatabase indexes using ArcPy v10.6, but I encounter an ERROR 000837 The workspace is not the correct workspace type.

Question: If arcpy.RebuildIndexes_management() works only with Enterprise geodatabases, then how to rebuild indexes in file geodatabases?

import arcpy, os

workspace = "C:/tmp/myDatabase.gdb"
fc = "MyFeatureClass"
arcpy.Exists(workspace)
>>>True
arcpy.Exists(os.path.join(workspace, fc))
>>>True

arcpy.RebuildIndexes_management(workspace, "NO_SYSTEM", [fc], "ALL")
>>>ERROR 000837 The workspace is not the correct workspace type

.


texts:
I am using ArcMap 10.6.1 and I am trying to export a map as a TIFF file. When I'm in layoutview the map looks ok, but when I export it, new layers show up. Or the scale of the layers change or something.

As you can see in the screenshot there are 7 citynames in the map. When I export it the layers change. It seems as the map scale is changing when I export it. I've checked the maximum and minimum scale of the layer and it is not what's causing the problem. I'm quite new to ArcMap, but this wasn't a problem a while ago.When I try to export old maps that worked fine before,I get the same problem. Could it be a function I've checked or something?


texts: Adding featureOver on map for multiple layers using carto.js and leaflet

texts:
How can I add featureOver on map for multiple layers(layers added dynamically) using carto.js and leaflet. This functionality have worked only for one layer. Here is my code:

layer.on(carto.layer.events.FEATURE_OVER, featureEvent => {
  //console.log("feature event", featureEvent);
  popup.setLatLng(featureEvent.latLng);
  if (!popup.isOpen()) {
    popup.setContent(featureEvent.data.sdlbl);
    popup.openOn(map);
  }
});
layer.on(carto.layer.events.FEATURE_OUT, featureEvent => {
  popup.removeFrom(map);
});


texts: Google Earth Engine: reproject images to the same projection

texts:
Below, I have clipped the same region from 2 google earth images. However, the two clips are at different projections (the crs argument in reproject determines the projection). The reason I am reprojecting is to decrease the resolution of the images, via the scale argument, so that they can be imported with sampleRectangle.  How would I clip the images at the same projection such that they overlap?
import ee
import matplotlib.pyplot as plt
from matplotlib.pyplot import imshow
import numpy as np

# Load images
pH = ee.Image("OpenLandMap/SOL/SOL_PH-H2O_USDA-4C1A2A_M/v02")
AAFC = ee.ImageCollection('AAFC/ACI')
AAFC = AAFC.filter(ee.Filter.date('2016-01-01', '2016-12-31')).first()

# Get geometry/polygon of interest
geometry = CAR.getInfo()['features'][15]['geometry']['coordinates']
region = ee.Geometry.Polygon(geometry)

# Clip AAFC and plot
AAFC_crs = AAFC.clip(region).projection().wkt()

AAFC_region = AAFC.clip(region)
AAFC_region = AAFC_region.resample('bilinear').reproject(crs=AAFC_crs, scale=1000)
AAFC_img = AAFC_region.sampleRectangle(defaultValue=0)
AAFC_img = np.array(AAFC_img.getInfo()['properties']['landcover'])
AAFC_img = AAFC_img.reshape(len(AAFC_img), -1)

plt.imshow(AAFC_img)
plt.colorbar()
plt.show()

# Clip pH and plot
pH_crs = pH.clip(region).projection().wkt()

pH_region = pH.clip(region)
pH_region = pH_region.resample('bilinear').reproject(crs=pH_crs, scale=1000)
pH_img = pH_region.sampleRectangle(defaultValue=0)
pH_img = np.array(pH_img.getInfo()['properties']['b30'])
pH_img = pH_img.reshape(len(pH_img), -1)

plt.imshow(pH_img)
plt.colorbar()
plt.show()


texts: Field Calculator by Expression Incorrect

texts: I've created a new field as a whole integer with a length of 10. I am attempting to use the field calculator to assign values to the newly created field. However, I am certain the numbers aren't being calculated correctly and I'm not sure why. At this point I'm just attempting to do 36*1000000000 and the result QGIS gives is 1640261632.


texts:
I'm currently running an application using Norwegian and Swedish basemaps, served by the respective government agencies via WMTS, but licensed under variants of Creative Commons. However, I'm not impressed by the latency and robustness of their solution. The Swedish map is slow in Norway and vice versa, and I can't even fathom how it looks over a trans-atlantic connection.

I'd like to copy the tiles of the relevant areas from the respective WMTS services and serving them myself over a CDN. Are there any industry standard tools for ripping the tiles of a WMTS or will I have to write my own script?


texts: Why are map scales not used instead of Zoom-levels in web mapping clients such as OSM or Google Maps?

texts: Google Maps and OpenstreetMap and many other web mapping products utilize the concept of zoom-levels to change the size of the map in a web client. Most people are however more familiar with the concept of map scales(i.e 1:1000, 1:50000 etc) when zooming in  or out of a map or considering measurements. When using zoom-levels the scales are not very intuitive (http://api.geoext.org/1.1/examples/zoom-chooser.html). Even desktop GIS software maintains the concept of scales when zooming in or out of digital maps. So why is the concept of map scale not used instead of Zoom-levels in web mapping clients such as OSM or Google Maps? This would enable the maps to be much more functional for administrative purposes for example in land use planing where the use of set scales like 1:10,000 are required when printing to paper.


texts: Dividing streets layer into 2 mile segments?

texts:
I have a streets layer which I want to divide into 2 mile segments. I have searched and read enough for days but there is no information on how to do this. The split command allows you to select and divide the poly-line which is manual. It is practically impossible for me to do it for a streets layer which has tens of thousands of poly-lines.

Is there any tool, script or combinations of tools I can use to achieve this?


texts: Location API works, but neither Google Maps, not Here maps can not find location

texts:
I used Holux M-241 as GPS receiver.
On Windows 8.1 I installed GPSDirect drivers, and now can see my location on PC.
I tested  that Location API can access to gps data (I used this Mozilla Developer Site to test)
However, no browser-based map application can not find my position. Why?


texts: Add bands from one dataset to another in Google Earth Engine

texts:
I want to add bands from one data set to another. I have processed two data sets (filterMOD09 and MODIS13). The resolution and time of their images are in one-to-one correspondence. What should I do?
var MOD = ee.ImageCollection("MODIS/006/MOD13Q1");
var MYD = ee.ImageCollection("MODIS/006/MYD13Q1");

var combinedImageCollection = MOD.merge(MYD);
var MODIS13 = combinedImageCollection
    .filterDate('2000-01-01', '2020-12-31')
    .sort("system:time_start")
    .select(['NDVI','EVI','sur_refl_b02'],['NDVI','EVI_Old','NIR'])
    .map(function(image){
      image = image.reproject('SR-ORG:6974',null,250);
      return image;
    });

var MOD09 = ee.ImageCollection("MODIS/006/MOD09A1");

var temporalCollection = function(collection, start, count, interval, units) {
  var sequence = ee.List.sequence(0, ee.Number(count).subtract(1));
  var originalStartDate = ee.Date(start);
  return ee.ImageCollection(sequence.map(function(i) {
    var startDate = originalStartDate.advance(ee.Number(interval).multiply(i), units);
    var endDate = originalStartDate.advance(
      ee.Number(interval).multiply(ee.Number(i).add(1)), units);
    return collection.filterDate(startDate).first();
  }));
};

var modis00 = temporalCollection(MOD09, '2000-02-18', 46, 16, 'day')
              .filterDate('2000-02-18','2000-12-31');

var modis01 = temporalCollection(MOD09, '2001-01-01', 46, 16, 'day')
              .filterDate('2001-01-01','2001-12-31');

var modis02 = temporalCollection(MOD09, '2002-01-01', 46, 16, 'day')
              .filterDate('2002-01-01','2002-06-27');

var combined_1 = ee.ImageCollection("MODIS/006/MOD09A1")
          .filterDate('2002-06-27','2020-12-31');
var combined_2 = modis00.merge(modis01).merge(modis02).merge(combined_1);

var filterMOD09 = combined_2
    .sort("system:time_start")
    .map(function(image){
      image = image.reproject('SR-ORG:6974',null,250);
      return image;
    })
    .select(['sur_refl_b06'],['SWIR']);


texts: ArcGIS 10 add-in: Bring dockable window to front

texts:
I have an ArcGIS 10 ArcMap add-in written in C# .NET 3.5 that implements an ESRI.ArcGIS.Desktop.AddIns.DockableWindow (which also inherits from UserControl) and an ESRI.ArcGIS.Desktop.AddIns.Tool that when clicked in the map updates the dockable window.

I'd like to bring the dockable window to the front of the Z-order in the Tool's OnMouseDown() method (when in undocked mode). Currently, if the user opens another dockable window, puts it on top of mine, and clicks with the tool, the window updates but it is not brought to the front. I am already calling IDockableWindow.Show(true) to ensure the window is visible after clicking with the tool. I have also tried UserControl.BringToFront() but it does not have any effect.

The best workaround I currently have is to call IDockableWindow.Show(false) followed by IDockableWindow.Show(true), however this not ideal as it is jarring to have the window disappear and reappear, as well as have to fully repaint which takes a substantial amount of time.

The built-in Identify window does not have this issue and is brought to the top each time the Identify tool is used, so clearly there is a way to do it.

Does anyone know of a better solution to this? Thanks!

Edit: Here is the code I ended up using to resolve this. Thanks Kirk and Petr!

public static void BringDockableWindowToFront(IDockableWindow dockableWindow, IntPtr dockableWindowControlHandle)
/// <summary>
/// Workaround for bringing a dockable window to the top of the Z order.
/// dockableWindowControlHandle is the Handle property of the UserControl implemented by the dockable window
/// </summary>
{
    IWindowPosition windowPos = dockableWindow as IWindowPosition;
    IntPtr parentHwnd = GetParent(dockableWindowControlHandle); // Get parent window handle
    switch (windowPos.State)
    {
        case esriWindowState.esriWSFloating:
            IntPtr grandParentHwnd = GetParent(parentHwnd); // Get grandparent window handle
            SetActiveWindow(grandParentHwnd);  // Activate grandparent window when in floating (undocked) mode
            break;
        //case esriWindowState.esriWSMaximize: // Mode not yet implemented in ArcGIS 10, check at 10.1
        //case esriWindowState.esriWSMinimize: // Mode not yet implemented in ArcGIS 10, check at 10.1
        case esriWindowState.esriWSNormal:
            SetActiveWindow(parentHwnd); // Activate parent window when in normal (docked) mode
            break;
    }
    SetFocus(dockableWindowControlHandle); // Set keyboard focus to the dockable window
}

// Retrieves a handle to the specified window's parent or owner.
[DllImport("user32.dll", SetLastError = true, CharSet = CharSet.Auto, ExactSpelling = true)]
public static extern IntPtr GetParent(IntPtr hWnd);

// Sets the keyboard focus to the specified window.
[DllImport("user32.dll", SetLastError = true, CharSet = CharSet.Auto, ExactSpelling = true)]
public static extern IntPtr SetFocus(IntPtr hWnd);

// Activates a window.
[DllImport("user32.dll", SetLastError = true, CharSet = CharSet.Auto, ExactSpelling = true)]
public static extern IntPtr SetActiveWindow(IntPtr hWnd);


texts: When I reconcile a version in Arc Map, if it has any conflicts it will stop the reconcile operation and open a window. The window shows the conflicts of that particular version with red dots. My task is to read those values and generate a report which contains only conflict fields.


texts:
I was working with a large data set that had data for roads, rivers, railways, etc. all together so I decided to break the different data out so I could label it differently on the map and check each on and off separately.  I selected all of one type, copied it to the clipboard, and pasted it as a new memory vector layer.  Everything was working fine yesterday, but when I opened the file this morning all of the data in the attribute table was gone in the vector memory layers I created.

I had toggled off edits and saved all my layers before closing yesterday.  Why did this happen? Should I have pasted as just a vector layer instead or should I have saved my selections as separate shapefiles instead?

Also, another issue I am experiencing is that I have several separate shapefiles with different address points.  I was copying and pasting points between shapefiles yesterday without issue (they are the same formatting in the attribute table btw) but when I tried to do that this morning I could paste, but not save the shapefiles.  I was getting the following error message:

  Errors: ERROR: 1 feature(s) not added.

  Provider errors:
        OGR error creating feature -2: Attempt to write non-point (MULTIPOINT) geometry to point shapefile.

Why was this working yesterday and not today? Did I somehow change the type of shapefile (i.e point, line polygon) without realizing it?

I am using QGIS 2.4.0-Chugiak.


texts: How to get the nearest line to a feature?

texts:
I have features on map, and i want to get the nearest poly-line to each one.
How to write spatial query to get the nearest one.

I tried to use ISpatialFilter but i couldn't reach anything.


texts: Where can the GDAL executables be found in a QGIS installation?

texts:
I run QGIS 1.8.0 Lisboa on a Windows 7 machine.  Yesterday I attempted to install the TauDEM plugin using these instructions: http://lists.osgeo.org/pipermail/qgis-developer/2012-June/020807.html

Since doing this the GDAL tools for QGIS do not work with error:

  The process failed to start. Either the invoked program is missing, or
  you may have insufficient permissions to invoke the program.

I understand this means they cannot find the GDAL executables.  However, neither can I.  Where actually are the GDAL executables in the QGIS folder?  I assume I then have to add this manually to the Environmental Variables or into the GDAL tool settings in QGIS - is this correct?

PS. I can run GDAL tools through the OSGEO4W Shell so GDAL itself is alive and well.


texts: I have rotated some text labels by -90 in my project.  I want them to stand vertically on an alignment bar.  If the map window with the alignment bar also has a rotation the text is no longer vertical or horizontal.  I am not sure where this text angle is coming from.  (Please see picture.)  The "3+200.00" is set horizontal and transposes between map windows correctly.  The "Dugout 4" is set to -90 and acquires an off-angle somewhere when it transposes to the map window.  Does anyone have any idea how to overcome this?


texts: custom plugin not showing - ImportError: No module named resources

texts:
I just develop a custom plugin in QGIS. It folder located here

  C:\Users(user).qgis2\python\Plugins\PluginTest

But it is not showing the plugin in QGIS, although a downloaded plugin VectorCompare showing correctly in plugin dropdown. I have no warning or error log.
I just saw an error in Plugin > Installed tab.

Couldn't load plugin TaskAutomation due to an error when calling its classFactory() method

>     ImportError: No module named resources
>     Traceback (most recent call last):
>       File "E:/QGISSE~1/apps/qgis/./python\qgis\utils.py", line 306, in startPlugin
>         plugins[packageName] = package.classFactory(iface)
>       File "C:/Users/FK/.qgis2/python/plugins\TaskAutomation\__init__.py", line 34, in
> classFactory
>         from .TLL_Task import TaskAutomation
>       File "E:/QGISSE~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
>         mod = _builtin_import(name, globals, locals, fromlist, level)
>       File "C:/Users/FK/.qgis2/python/plugins\TaskAutomation\TLL_Task.py", line 26, in
>         import resources
>       File "E:/QGISSE~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
>         mod = _builtin_import(name, globals, locals, fromlist, level)
>     ImportError: No module named resources


texts: Transforming objects to feature collections in Google Earth Engine

texts:
I'm using GEE to export a panel dataset with population estimates and forest loss in Rwanda 2000-2011. My code is below:

//Rwanda admin shapefile
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

var Rwanda = countries.filter(ee.Filter.eq('country_co', 'RW'));
Map.addLayer(Rwanda);

//Hansen Global Forest Change data
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015')

//Population data
var Allpop = ee.ImageCollection("WorldPop/GP/100m/pop")
  .filterDate('2000', '2012');

var pop = Allpop.filterMetadata('country', 'equals', 'RWA');

//////////////////////////////////////////////////////////////////
//Annual forest loss in Rwanda using the GFC dataset
//////////////////////////////////////////////////////////////////

var lossImage = gfc2014.select(['loss']);
var lossAreaImage = lossImage.multiply(ee.Image.pixelArea());

var lossYear = gfc2014.select(['lossyear']);
var lossByYear = lossAreaImage.addBands(lossYear).reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1
    }),
  geometry: Rwanda,
  scale: 30,
  maxPixels: 1e9
});

var statsFormatted = ee.List(lossByYear.get('groups'))
  .map(function(el) {
    var d = ee.Dictionary(el);
    return [ee.Number(d.get('group')).format("20%02d"), d.get('sum')];
  });
var statsDictionary = ee.Dictionary(statsFormatted.flatten());
print("Annual forest loss", statsDictionary);

//////////////////////////////////////////////////////////////////
//Annual population in Rwanda using the WorldPop data
//////////////////////////////////////////////////////////////////

function newCollectionToImage(collection){
var stack = ee.Image(collection.iterate(function(img, prev) {
return ee.Image(prev).addBands(img);
}, ee.Image(1)));

stack = stack.select(ee.List.sequence(1, stack.bandNames().size().subtract(1)));
return stack;
}

var tch = pop.map(function(img){
var y = ee.String(ee.Date(img.get('system:time_start')).get("year"));
return img.set("year",y).rename(y);
});

var tchImg = newCollectionToImage(tch).clip(Rwanda);
var pop_Rwanda=tchImg;

var Rwanda_annualpop = pop_Rwanda.reduceRegions(Rwanda, ee.Reducer.sum());

print("population", Rwanda_annualpop);

The variable statsDictionary (annual forest loss) is an object, and there is very little documentation online about how to transform these objects into something that can be exported.

I would really like to export one dataframe that includes, for each year, the population estimate and the forest loss. Is there a way to transform the object into a feature collection to facilitate this?


texts:
I am trying to model property line setbacks. In this case: 25 ft from front and back of property, 10 ft from one side and 8 ft from the other. Tried averaging the total set back distance 68/4 = 17 ft from each side. However this does not produce the same result as the specified setbacks. Typically off by 100-200 sq. ft.

Are there any tools similar to buffer that will model the specific setbacks. I could identify the centroids and somehow specify from there, or use the street centerlines.


texts: Google Earth imagery not displayed using SAS.Planet:

texts:
I have some trouble using SAS.Planet to download some georeferenced Google Earth images. If I try to open (display) "Earth(GeoCacher direct)" (or basically any other Google Earth option) there are no tiles building up. I just get a grey background. Also there is no Errorcode shown. I updated the Mapversions without any changes. SAS.Planet Version: SAS.Planet.Nightly.201020.10106
Is there a way to get the Google Earth Imagery working in SAS.Planet?


texts:
I am trying to execute this code,

import osgeo.ogr

def findPoints(geometry, results):
    for i in range(geometry.GetPointCount()):
        x,y,z = geometry.GetPoint(i)
    if results['north'] == None or results['north'][1] < y:
        results['north'] = (x,y)
    if results['south'] == None or results['south'][1] > y:
        results['south'] = (x,y)
    for i in range(geometry.GetGeometryCount()):
        findPoints(geometry.GetGeometryRef(i), results)

shapefile = osgeo.ogr.Open("../../Data/tl_2009_us_state/tl_2009_us_state.shp")
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(53)
geometry = feature.GetGeometryRef()

results = {'north' : None,
           'south' : None}

findPoints(geometry, results)

and I am constantly receiving this error,

Traceback (most recent call last):
  File "identify_northsouth_point.py", line 22, in <module>
    findPoints(geometry, results)
  File "identify_northsouth_point.py", line 8, in findPoints
    results['north'] = (x,y)
UnboundLocalError: local variable 'x' referenced before assignment

I have tried global and nonlocal, but it does not work. Since I am not gettin any input from outside the function, so I anyways would not require global or nonlocal.


texts: Publishing tif with elevation data to geoserver

texts:
I Have been trying to publish a DEM built from SRTM data (both a vrt of the .hgt files and a GeoTiff converted from the vrt).
I Keep getting Cannot enable, no attribute of type Number found on both elevation and time dimensions no matter what data type i save them in.
Is there a a specific data type which GeoServer will recognize as Number?
If so what is the conversion method that will work on best on the tif file with GDAL or QGIS?


texts: Updating tags on user's portal item via API

texts:
I'm using the extract data widget to create extract files.  It writes these files to the logged in portal user's content.
I want to tag the items that were just created but I can't find anywhere in the api that lets me do that.
Example of what I want to do:
var myUrl = 'somevalue';  // I know this value at runtime, so not a problem.
var myItemId = 'someOtherValue'; // Ditto.

var thisPortal = new arcgisPortal.Portal(myurl);
thisPortal.signIn().then(function (loggedInUser) {
    loggedInUser.getItem(myItemId).then(function(theItem) {
         // I want to add a tag to the item so I know it's ok to automatically toss it a day later.
         theItem.tags[theItem.tags.length] = 'Delete in One Day';
        // I want to write the item back to the portal.
        // And I cannot find any documentation on this step!
    })
});

I don't want to have to trust the users to clean up after themselves.   I have a concern that they will clog up the disk with extract files that they will never look at again and never throw away again.
I can ask the users to mark whether the extract should be held temporarily or permanently before the extract files are created (and default to temporarily!)  If I then tag the extract files accordingly, I can have an automated program run every evening and toss out the expired ones.


texts: open source command line tools on Linux or source code for converting tiff image file to BSQ file

texts: Let me know homepages related to open source command line tools on Linux or source codes for converting tiff image file to BSQ file.


texts:
I have created a new control from an inheritance:

L.Control.extend({
            options: {
                position: position,

                descriptionHeader: 'This is a description',
                descriptionText: 'Label',

                counterHeader: 'This is a description',
                counterText: 'Label',

                className: 'leaflet-control-custom',
                labelClassName: 'button-label'
            },
            setDescriptionHeader(text: string){
                console.log("here");
                this.options.descriptionHeader = text;
            },
            setDescriptionText(text: string){
                this.options.descriptiontext = text;
            },
            setDescriptionCounterHeader(text: string){
                this.options.counterHeader = text;
            },
            setDescriptionCounterText(text: string){
                this.options.counterText = text;
            },
            onAdd: function (map) {
                let mainContainer: any = L.DomUtil.create('div', 'description-control');
                let container: any = L.DomUtil.create('div', 'description-options');
                let counterContainer: any = L.DomUtil.create('div', 'description-options-2');

                let header: HTMLElement = document.createElement("h2");
                header.className = 'description-title';
                header.innerText = this.options.descriptionHeader;
                container.appendChild(header);

                let text: HTMLElement = document.createElement("span");
                text.className = 'description-text';
                text.innerText = this.options.descriptionText;
                container.appendChild(text);

                let counter: HTMLElement = document.createElement("h2");
                counter.className = this.options.labelClassName;
                counter.innerText = this.options.counterHeader;
                counterContainer.appendChild(counter);

                let counterText: HTMLElement = document.createElement("span");
                counterText.className = this.options.labelClassName;
                counterText.innerText = this.options.counterText;
                counterContainer.appendChild(counterText);

                mainContainer.appendChild(counterContainer);
                mainContainer.appendChild(container);

                return mainContainer;
            }
        });

I want this control content to be updated automatically, when an external variable change or when an action is triggered:

I want to be able to pass some variables to the control to update its text.
I've tried several solution like setting setters for variables, but its not working.

The control looks like the following:

How do I update the content of my control?


texts: How to transfer a WKT string from OpenLayers to Geoserver or PostGIS?

texts: Supposing I have created a Vector Layer in OpenLayers and then transformed it to WKT or GeoJson, how can I export it from OpenLayers, and possible import it to PostGIS or GeoServer?


texts:
I have been struggling with a CSV file I need to import into QGIS 2.2.0 (on Windows 7).

The coordinates are in this format:

LATITUDE  LONGITUDE
424605    861204
395222    861206
303051    861219
32155     861220
395541    861229
425119    861242

Latitude and Longitude are in units of degrees, minutes and seconds according to the person who provided the data.

I tried different formats (placing commas, dots, ',", etc at different places) to import this dataset into QGis, checking the "DMS coordinates" box in the "Create a layer form Delimited Text File". I always had a Delimited text file error message.

What is the format expected by QGis for Latitude/Longitude in DMS?


texts: Paths not matching on map, while they do in Illustrator

texts:
I'm quite a newcomer when it comes to map data.
I have an Illustrator map, drawn by myself, with different polygons (regions, country, ...).
In Illustrator I know 100% that the region and country line are perfectly in the same position.

I then used MAPublisher to export these to KML files, to then convert them into GeoJSON files.
I use these GeoJSON files together with Leaflet to generate polygons in the shape of what I have drawn in Illustrator.
But as you can see in the following picture, when I zoom in, the lines of regions and country are not in the same position.

Any ideas why this is the case?


texts:
I'm not terribly familiar with ArcGIS in general, but I have a geodatabase (I think) containing a 3m DEM in ArcGis 9.3.

I would like to have this 3m DEM in a GeoTIFF file so I can do some processing with the TauDEM package.

I've tried right-clicking the 3m DEM layer and going to Export->Tiff, but it asks me what size (DPI, in essence) I want the exported TIFF to be.

If I put in a low DPI, the exported file size is too small to be the actual DEM and comparisons between the GeoTIFF image and my view on ArcGIS show that the TIFF is smoothed. If I put in a high DPI I can't open the TIFF, but I'm pretty sure the result wouldn't be what I wanted, based on the low-DPI test.

How do I export a GeoTIFF containing exactly the data in the 3m DEM layer without having to worry about resolution?

Thanks!


texts: I want to extract the Z coordinates from the line nodes. The data is there in the vector editor.


texts: While connecting MS SQL getting error "Data source name not found and no default driver specified QODBC3: Unable to connect"

texts:
I have used the following code for connecting the Microsoft SQL Server's database. but at the time of opening database connection it is throwing an error.

Here is my code,

 from PyQt4.QtGui import *
 from PyQt4.QtSql import *
 from PyQt4.QtCore import *
 import sys

 db = QSqlDatabase.addDatabase("QODBC")
 db.setHostName("DB\DBCENTER2008R2");
 db.setDatabaseName("AU_Test");
 db.setUserName("sa");
 db.setPassword("sa12");
 if (db.open()==False):
     QMessageBox.critical(None, "Database Error",db.lastError().text())


texts:
The new R-bridge, unveiled at ESRI UC conference (2015), has a dependency on ArcGIS 10.3.1.

Has anybody tested it with older versions of ArcGIS Desktop? Specificially, I need to know why it depends on 10.3.1. I need to run it on 10.2.


texts: Use GDAL to resample and stack multidate image spectral bands

texts:
I have a Sentinel 2 multidate image data set and I'm trying to resample all of them at 10m and then stack them to apply a cloud mask, in Python language.
Being new to programming, I have a hard time getting started. By cons I have some ideas like:
  - recover the coordinates (x, y) of one of the 4 spectral bands at 10m (B2, B3, B4 or B8) and apply them on all the other bands of the multidate images.
The fact is that the ultimate goal is to create an atlas, so I also asked the question of whether to create a stack for each image in a function, but it looks very complicated. ..


texts: Project data into WGS from rotated raster (CMIP5 file ) in R

texts:
I have this raster (less than 10 mb). It is a historical simulation model (SMHI) created by CMIP5 clipped for the city of Amsterdam.

dropbox <- here is the file

What I want to do is to convert this rotated lat/lon to a WGS projection and then re-write the newly obtained raster using ncdf4 or Raster package in R (preferably ncdf4 since I have to do this with many files)

How to do this exactly?

I was trying something like for a single layer:

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

rs <- brick("SMHI_MPI-M-MPI-ESM-LR_historical_r1i1p1_RCA4_amsterdam_1970_2005.nc", varname = 'pr')

rp <- rasterToPoints(rs[[1]])
rp[,3] <-  rp[,3] * 86400 # the precipitation is measured in kg m's (convert here to daily mm)

# change the proj
r2 <-project(rp[,1:2], paste("+proj=ob_tran +o_proj=longlat", "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"), inv=TRUE, use_ob_tran=TRUE)

# add the projected coordinates to the matrix; convert it to df
rp[,1:2] <-r2
rp <-as.data.frame(rp)

# convert it to spatial
coordinates(rp)<-~x+y

# plot it
spplot(rp, scales = list(draw = T))


texts: Extracting pixel values in satellite images using ArcGIS for Desktop?

texts:
I'm trying to extract the pixel values(intensity values) of Landsat-8 satellite images using ArcGIS, I have an image and a crop data layer(CDL).

Then I only kept the data for one crop, and got the following result

Now I want to know pixel intensity for these selected areas on the original satellite image.

I combined images together and  use "Pixel Inspector" to see the values.

But I want to have all the values, so I can then analyse them and find some specific values like max or min.

What I need is the pixel intensity in the satellite image but only in those areas which are highlighted with green in "green and white Image" , so I need to find the location of those green areas then go to satellite image and find the intensity value in the corresponding location.

How can I do that?


texts: Featureclass attribute tables corrupted, but still exportable?

texts:
I have a model that takes a geodatabase of rasters and a featureclass, iterates through the rasters using the FC to calculate zonal stats on each, then creating a field and storing the values in the FC.    It actually works as intended, though it does have an unhandled python error when I compute a field and try to divide null by a float.  My question is about a really weird problem I encountered when I let it run overnight.  I queued a couple dozen iterations of this model and left them to run overnight.  Sometime in the night, my VPN disconnected from my work network, which caused ArcGIS to lose its connection to the license server and halt operations.

The next morning, I saw this, reconnected the VPN, and it seemed to sluggishly resume the model runs it had not completed before the disconnect, ultimately completing without errors.  However, upon looking at the results, everything but the first few were heavily corrupted.  My guess is these are the ones that completed before the disconnect happened.  For the corrupted ones their geometry is all still intact (draws properly), but the attribute tables are trashed.  They either will not open:  or have completely mangled values:
For reference, 'Name' should contain a place name and 'prod' should be a number somewhere in the 5-9 figure range.

When I realized what happened, I started replacing the corrupted feature classes to start the whole process over. During this process I realized that if I export features from one of the corrupted files, I get the feature class in the state that it was in before the model run (uncorrupted, but without all the zonal stats loaded into new fields).  This tells me that the data is still there.  Each of the model runs gave output saying that it was attaching the data to the FC, so presumably, that data should be there as well?

Does anyone understand what happened on a technical level?  I assume this has something to do with having the files open for writing and not saving changes properly due to the unexpected halt, but I would have thought ESRI would have planned for that and exited cleanly.  I'd like to try to recover the data, but I'm not sure where to start as I don't really understand exactly why they got corrupted.


texts: Split polygon with polygon in different tables / layers

texts:
I have a large number (1 million+) multi point polygons (red lines, table_a), in which are around 20 million other multi point polygons (black blocks table_b). They are in separate tables. What i want to do is measure the sum of the area of the black blocks within the red line.

My workflow is to first cut the black blocks along the red line then use st_intersection to select them and sum using st_area.

Something like

SET table_a
 UPDATE bldg_footprint = st_area (cut_buildings)
 FROM table_b
 WHERE ST_Intersects( admin_bound table_a , buildings table_b )

I can see from this answer here. That I can either the blocks directly using the polygon or first extract the lines and then cut using lines.

My question is wether there's a function in PostGIS that is more efficient than that approach?

I'd be open to a QGIS solution but had assumed that given the size of data PostGIS will be faster.


texts: What LiDAR processing tools are available in Python?

texts:
I have been using FUSION and the command line FUSION Lidar Toolkit (LTK) to process LiDAR data. A broad Google search ("Lidar Python") yielded libLAS and pyLAS as Python LiDAR libraries, however, these appear to provide only read and write access to LAS data.

I am particularly interested in creating intensity and density images in addition to canopy surface models from point clouds. Is there a generally accepted set of tools in Python that can accomplish the same sort of tasks FUSION LTK is capable of?


texts: Create raster surface from points with barrier

texts:
I have a point layer, that represents an elevation in that given point. This I need to create a raster/surface/DEM from, preferably by interpolating by Nearest Neighbour method or some way Liniar between the points. However I only need it inside my barrier which I have as a polygon and polyline layer. I'm trying to use Nearest/Natural Neighbour method with the polygon/polyline layer as processing extent, but this doesn't give the "cut" raster result I'm looking for.

I've tried converting my buffer into a raster and add that under Processing Extent Snap Raster, without getting the desired results.

I have ArcGIS Advanced 10.5.1 with Spatial and 3d Analyst.


texts: Make feature layer based on date, problem in python

texts:
I have a script that runs, that takes the current date and puts into a variable, which is then used to calculate the date field for the Feature class. Further in the script, the updated data is loaded into the live SDE feature class.

The old data get deleted with the following task.

1 - Make Feature layer, where the date field is not equal to the date variable.

2 - Delete these features.

My problem is that the to create feature layer is empty as i can't get the syntax right.

I have tried the following

gp.MakeFeatureLayer(myOutputConnexion + os.sep + "SY_PLN_MON_INFO.LLPG_POINTS","OLD_LLPG_POINTS","GIS_TIMESTAMP <> TO_DATE('" + ExactScriptTime + "00:00:00' ,'YYYY/MM/DD HH24:MI:SS')")

and

    gp.MakeFeatureLayer(myOutputConnexion + os.sep + "SY_PLN_MON_INFO.LLPG_POINTS","OLD_LLPG_POINTS","GIS_TIMESTAMP <> TO_DATE(' + ExactScriptTime ','YYYY/MM/DD HH24:MI:SS')")

The ExactScriptTime variable value is 21/08/2012 13:00:07 (which gets created in at the start of the script).

I know that the feature layer is empty because i have exported the feature layer into an FC, and it is empty.

Has anyone got any ideas on where I am going wrong.

Thanks


texts: QGIS: Legend with more columns in one row

texts:
I have two layers (water_supply_design, water_supply_condition). I'd like to create new map with Map Composer in QGIS. In this map, I'd like to have a legend with one row and two columns (Water supply, design, state). Is it possible to generate from this two layers a legend in that form?

Is this possible to create in the map composer?


texts: ArcGIS Calculate Field, Python Shape Expression

texts:
Using ArcGIS 10.1 and Python 2.7-
I am trying to write an expression that calculates the mid point of a line.  I've tried many different expressions, but I continue to get syntax errors.

#An example of one of my many failed attempts
arcpy.CalculateField_management(StreetCenterlines, "Latitude", "!SHAPE@X!")
arcpy.CalculateField_management(StreetCenterlines, "Longitude", "!SHAPE@Y!")

Can anyone lend a hand at how my expression should read in order to accomplish this?

The error I receive is: ERROR 999999: Error executing function. Syntax error Failed to execute (CalculateField).


texts: ERA5 grib file: how to know what each band means?

texts:
I downloaded monthly averaged data for a given area, 2 variables and 3 years. I read the grib file on R and I see there are 72 bands (12*2*3). Great. But how do I know what each one is? I would expect some kind of legend (like this). But nothing in the Spatial Grid Data Frame seems recognizable.

Any idea?

The file can be dowloaded here.


texts: QGIS Server [Map request failed error: Host _ not found]

texts:
I am trying to set up a QGIS Server to publish the WMS layer using a docker container.
I followed the following section in the QGIS Manual https://docs.qgis.org/3.16/en/docs/server_manual/containerized_deployment.html#docker-stacks
Based on this I was able to see this page successfully using the following url and it looked like that everything had been set up correctly. http://localhost:8080/qgis-server/?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

I was also able to add the WMS connection and see layers in QGIS Browser like this

But when I try to load this layer in the canvas I get the following error

I tried publishing the GeoJson Polygon layer by adding it to the source QGIS project but it was giving the same error as well. What am I doing wrong?


texts: PyQGIS: New vector created inside group instead of TOC

texts:
def newPolygonVector(self):
    vlt = QgsVectorLayer("Polygon?crs=EPSG:32647", "", "memory")
    vlt.dataProvider().addAttributes([QgsField("id", QVariant.Int)])
    vlt.updateFields()
    self.MakeFile()
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.driverName = "ESRI Shapefile"
    QgsVectorFileWriter.writeAsVectorFormatV2(vlt, self.mypath + f'/FieldBorderPolygon{self.index}.shp',
                                              project.transformContext(), options)
    self.vlb = self.iface.addVectorLayer(self.mypath + f'/FieldBorderPolygon{self.index}.shp', "", "ogr")
    Symbol = QgsFillSymbol.createSimple(
        {'color_border': '#3232fa', 'width_border': '1', 'style': 'no'})
    self.vlb.renderer().setSymbol(Symbol)
    self.vlb.startEditing()
    self.vlb.triggerRepaint()

def newLineVector(self):
    vlt = QgsVectorLayer("LineString?crs=EPSG:32647", "", "memory")
    vlt.dataProvider().addAttributes([QgsField("id"
                                               , QVariant.Int)])
    vlt.updateFields()
    self.MakeFile()
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.driverName = "ESRI Shapefile"
    QgsVectorFileWriter.writeAsVectorFormatV2(vlt, self.mypath + f'/RajahLine{self.index}.shp',
                                              project.transformContext(), options)
    self.vlr = self.iface.addVectorLayer(self.mypath + f'/RajahLine{self.index}.shp', "", "ogr")
    self.vlr.renderer().symbol().setWidth(0.6)
    self.vlr.renderer().symbol().setColor(QColor.fromRgb(255, 0, 0))

def run(self):
    """Run method that performs all the real work"""

    # Create the dialog with elements (after translation) and keep reference
    # Only create GUI ONCE in callback, so that it will only load when the plugin is started
    if self.first_start == True:
        self.first_start = False

    self.newLineVector()
    self.newPolygonVector()

    #Grouping
    root = project.layerTreeRoot()
    group = root.insertGroup(0, f"Prelining {self.index}")
    group.addLayer(self.vlb)
    group.addLayer(self.vlr)
    root.removeLayer(self.vlb)
    root.removeLayer(self.vlr)

These are some of my code snippet for one of the processes in a plugin and the plugin is iterable.

Whenever it runs the first time, it gets the result that I want, as it creates a 2 new vector outside of a group in TOC

But when it comes to the second time I run the plugin, it will create the first 2 vector inside the 1st group and then will add it into 2nd group having duplicated layer in Group1 and Group2

Results are shown below

1.1) This is the first time running the plugin and everything lays greatly inside the group:

1.2) When I click on the iterable plugin again, it will create the vector inside the previous group instead of outside the group:

1.3) Making me having 2 duplicated layers:


texts: Calculating NDVI for Sentinel 2A, graph builder (SNAP)

texts:
While I am trying to calculate NDVI in graph builder, I get this error:

Can not convert "(B8-B4)/(B8+B4)" to double.

Target band type is: Float32 (it comes up automatically when I load the image.


texts:
I want to join the results of sDNA analysis to urban blocks and I have a list of questions
1- Do I need to weight sDNA values of the street segments (i.e. the links) before spatial join to urban blocks layer?
2- if so, what is the formula for the weighting variable?
3- How could I aggregate the sDNA values of surrounding segments of each block?
4- and lastely, could anybody in simple language illustrate the difference between betweenness and two phase betweenness and when should I depend on the two phase betweennness?


texts: How to find similar GPS tracks considering timestamps?

texts:
There are around 300 users tracks per day with about 4k points in each with 10s time interval. I'd like to find similar tracks considering location and time.
The tracks are similar if it can be assumed that two users where moving together (close locations during 30s timeframe) for the most of the track.

I have already seen similar rather old questions:
How to find if two gps tracks were from the same trip? and
Matching GPS tracks
but none of them answers what algorithm can be used find a time-based solution.

Maybe there are some new ideas after 5 years?


texts: Only uploading group layers instead of single simple layer when using ArcPy

texts:
I use the getWebLayerSharingDraft to upload layers online with Python. However, whenever I upload it, the result is always a group layer with a single layer inside. I wonder if it is possible to use Python to upload as a simple layer instead. My code:

arcpy.SignInToPortal('https://realtermenergy.maps.arcgis.com', username, password)
gis = GIS('https://mysite.maps.arcgis.com', username, password)
outdir = r"C:\Folder"
sddraft_filename = "name.sddraft"
sddraft_output_filename = os.path.join(outdir, sddraft_filename)

# Create FeatureSharingDraft and set service properties
sharing_draft = aprxMap.getWebLayerSharingDraft("HOSTING_SERVER", "FEATURE", GroupLayer)
sharing_draft.summary = "This is a layer."
sharing_draft.tags = "tags"

# Create Service Definition Draft file
sharing_draft.exportToSDDraft(sddraft_output_filename)

# Stage Service
sd_filename = "layername.sd"
sd_output_filename = os.path.join(outdir, sd_filename)
arcpy.StageService_server(sddraft_output_filename, sd_output_filename)

# Share layer to portal
arcpy.UploadServiceDefinition_server(sd_output_filename, "My Hosted Services")

The funny thing is that the Arcgis Online page is identical to layers upload without using python (using the sharing as a weblayer tool on arcgis PRO), yet the layer using python is a group layer, and using the Pro tool is a single one.

Whenever I go to the catalog on arcgis PRO, and search for the layer and bring to the map, it brings as a group layer named GroupLayer, with the layername Layer inside.

How I use a python code and upload a grouped layer?


texts: How to convert GeoJSON points to H3 hexagons in Python

texts:
I am trying to take a geojson file of points and convert it to H3 hexagons in my Python script. In JavaScript, there is a separate package called geojson2h3, but the Python h3 module has a method called `.geo_to_h3' that I am trying to utilize. It seems to take 3 arguments: a lat(X), a long(Y), and a resolution level. I have pulled in my GeoJSON file and have a JSON object like this:
{"crs": {"properties": {"name": "EPSG:4326"}, "type": "name"},
     "features": [
       {"geometry": {"coordinates": [-75.330278, 39.643889], "type": "Point"},
        "id": 0, "properties":
           {"FID": 0, "Gauge_Name": "Salem River At Woodstown Nj", "Gauge_No": "01482500",
            "State": "NJ", "X": 39.643889, "Y": -75.330278}, "type": "Feature"},
       {"geometry": {"coordinates": [-75.255556, 39.4725], "type": "Point"},
        ....

I can pull out coordinates of the individual features, but its messy and doesn't seem efficient:
import h3
import geojson
import os
import pandas
import numpy as np
import geopandas as gpd

dc_json_dir = r"C:\Projects\H3_Uber_Hexagon\data"
usgs_gauges = "riv_gauge_points.json"
gauge_path = os.path.join(dc_json_dir, usgs_gauges)

with open(gauge_path) as f:
    gj = geojson.load(f)
features = gj['features'][0]

x = features.geometry.coordinates[0]
y = features.geometry.coordinates[1]
level = 7

hexagon_test = h3.geo_to_h3(x, y, level)
print(hexagon_test)
>>>'87f05c0c0ffffff'

So, it returns the H3 spatial index, but I still need to add it to the GeoJSON object or an object that I can render as a hexagon later. I can also put the geojson into a dataframe and maybe add the data from the X, Y columns:
def hex_indexer(row):
  x = row.X
  y = row.Y
  level = 10
  hex_val = h3.geo_to_h3(x, y, level)
  return hex_val

gauge_df = gpd.read_file(gauge_path)
gauge_df['H3_Index'] = gauge_df.apply(hex_indexer, axis=1)
gauge_df.head()

    FID State  Gauge_Name      X          Y         Gauge_No    geometry                    H3_Index
  0   0 NJ     Salem River  39.643889  -75.330278   01482500    POINT(-75.33028 39.64389)   8a2aaca5132ffff
  1   1 NJ  Cohansey River  39.472500   -75.255556  01412800    POINT (-75.25556 39.47250)  8a2aadc34847fff

While this seems to work to create the h3 index, something doesn't seem to be right. Is there a better way to do this and get the data in a format that can be rendered using Python?


texts: FlexViewer 2.1 - Add scroll for widgetgroup

texts:
I have a site that is using the ArcGIS Viewer for Flex version 2.1.  In the header I have a widget group that contains 25 items.  The list is so long that it goes beyond the bottom of the browser window and I can't see the last 5-6 items.  Is it possible to add scroll capability to the widgetgroup so I can scroll through the items?

Thanks,
Tim


texts: Converting a FeatureCollection to an ImageCollection

texts:
I am trying to convert a FeatureCollection (a bit heavy one though) to an Imagecollection.
It seems that the job has been done successfully, however, I can't take any further steps with the Imagecollection such as clipping or reducing, and I cannot address the issue,
Ang thoughts on this?
Here is the link to the code:
https://code.earthengine.google.com/31de43e96b218cac08248e9177800d7e


texts: Network Analyst With Centroid Connectors

texts:
I have a fairly straight forward routing problem that I cannot seem to solve with the various Network Analyst tools available to me. We have categorized our street network according to the level of traffic stress one experiences while biking (colored lines in the image below). The goal of this analysis is to calculate the network distance between the top trip producing census blocks (green points in the image below) and top trip attracting census blocks (red point in the image below). This would be done twice, once for all levels of traffic stress and once for only minimal levels of traffic stress so a comparison can be made.

In order to "connect" the census block centroids to the network, I created connectors (the gray lines emanating from the points) and integrated them into our network. The number of connectors was based on how many network segment midpoints fell along the border of the census block. I have weighted these connects with 0 cost which was a decision made by our research team. The problem I am having, is that the route algorithm obviously searches for the route that costs the least (black line below) causing the route to traverse multiple centroid connectors and avoiding the network all together. I have introduced turn restrictions (see the red lines adjacent to the origin of the trip below) that restricts movement through two connectors. This appears to work, but would mean I would have to build turn lines for each pair of connectors in my network!

Does anyone know of an alternative solution to this problem? If there was a way to tell Network Analyst to never pass through a centroid or to only traverse one connector for any given route solution. I cannot think of a way to do this. Any ideas?


texts:
I am trying to look at how each input variables contribute to my Random Forest classifier. There is an explain() tools in ee.Classifier but the portal always shows:
 Dictionary (Error)
 Earth Engine capacity exceeded.

It is a dictionary, is there any way to export the explanation results so I can check it externally like in Excel?
var trainedClassifier = ee.Classifier.smileRandomForest(10).train({
    features: training,
    classProperty: 'LCZ',
    inputProperties: bands
    });

    var classified = composite.classify(trainedClassifier);

    print('explain', trainedClassifier.explain());

I don't know what code to write for exporting then.


texts: What causes a Crashed dumped error when adding labels to large file?

texts:
I created a world grid (rectangles) using mmqgis and joined a CSV file to the attribute table in order to label the grid.

Everything runs fine but when I try to add labels to the grid I get a minidump and QGIS shuts down.  Smaller sample sizes of the same data work fine.

Am I trying to do too much by adding 90,000+ labels?


texts: Calculating distances between n geographical points?

texts:
I would like to calculate the distances between N (an arbitrary number, but usually <50-100) geographical points (e.g. distance between points a, b, c, d, e, that is to say between a and b, a and c, a and d and so forth) in km (via latitude and longitude) and to draw a 2D bitmap map (PNG, JPG) of it all with the straight-line distances (in km) shown near the lines that represent them.

I'd prefer an online tool, but a FOSS GIS package is also an option. I'm talking about any geographical points on planet Earth. The map should be either topographic, satellite imagery, or combined (like Google Maps'/Google Earth's), either 2D or 3D. The inclusion of basic labels is good to have (country names, seas, oceans, rivers, lakes, dams, mountains, settlements, etc.).

I've found this wonderful tool: http://www.gpsvisualizer.com/calculators, but I'm not sure what projection it uses and if its calculated distances can be deemed rather (approximately) correct.

I know that Google Maps and Google Earth offer measurements between an N number of geographical points, but I haven't discovered how to set the points themselves in terms of latitude and longitude, you have to drag dots to your desired locations and I'm not sure how to see their longitude and latitude.

Can I do that using http://www.OpenStreetMap.org?

The topographic or satellite imagery level of detail of the map is not expected to be very high.


texts:
Getting the following error message when running r.viewshed

Processing algorithm…
Algorithm 'r.viewshed' starting…
Input parameters:
{ '-b' : True, '-c' : True, '-e' : False, '-r' : False, 'GRASS_RASTER_FORMAT_META' : '', 'GRASS_RASTER_FORMAT_OPT' : '', 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'coordinates' : '574583.6879638762,4750016.0627968935 [EPSG:26710]', 'input' : 'C:/Users/seansuzedielois/Documents/GIS/Mastering QGIS/Exercises/elev_clip.tif', 'max_distance' : 32000, 'memory' : 500, 'observer_elevation' : 20, 'output' : 'C:/Users/seansuzedielois/AppData/Local/Temp/processing_5a9b7eb19ebd43a9abc1f2f18e08d391/c11c8acdc82b4f71b2eb4c4ab3a15496/output.tif', 'refraction_coeff' : 0.14286, 'target_elevation' : 0 }

g.region n=0.0 s=0.0 e=0.0 w=0.0 res=100.0
r.viewshed input=rast_5d05e88ecd7383 coordinates=574583.6879638762,4750016.0627968935 observer_elevation=20 target_elevation=0 max_distance=32000 refraction_coeff=0.14286 memory=500 -c -b output=output18d498a864564b1ebfcce3671ea42038 --overwrite
g.region raster=output18d498a864564b1ebfcce3671ea42038
r.out.gdal -t -m input="output18d498a864564b1ebfcce3671ea42038" output="C:\Users\seansuzedielois\AppData\Local\Temp\processing_5a9b7eb19ebd43a9abc1f2f18e08d391\c11c8acdc82b4f71b2eb4c4ab3a15496\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite
Cleaning up temporary files...
Starting GRASS GIS...
WARNING: Concurrent mapset locking is not supported on Windows
Executing <C:\Users\seansuzedielois\AppData\Local\Temp\processing_5a9b7eb19ebd43a9abc1f2f18e08d391\grassdata\grass_batch_job.cmd> ...
C:\PROGRA~1\QGIS3~1.4\bin>chcp 1252 1>NUL
C:\PROGRA~1\QGIS3~1.4\bin>g.region n=0.0 s=0.0 e=0.0 w=0.0 res=100.0
ERROR: North must be larger than South
C:\PROGRA~1\QGIS3~1.4\bin>r.viewshed input=rast_5d05e88ecd7383 coordinates=574583.6879638762,4750016.0627968935 observer_elevation=20 target_elevation=0 max_distance=32000 refraction_coeff=0.14286 memory=500 -c -b output=output18d498a864564b1ebfcce3671ea42038 --overwrite
WARNING: East-west resolution does not equal north-south resolution. The viewshed computation assumes the cells are square, so in this case this may result in innacuracies.
Computing events...
ERROR: Raster map [rast_5d05e88ecd7383] not found
C:\PROGRA~1\QGIS3~1.4\bin>g.region raster=output18d498a864564b1ebfcce3671ea42038
ERROR: Raster map <output18d498a864564b1ebfcce3671ea42038> not found
C:\PROGRA~1\QGIS3~1.4\bin>r.out.gdal -t -m input="output18d498a864564b1ebfcce3671ea42038" output="C:\Users\seansuzedielois\AppData\Local\Temp\processing_5a9b7eb19ebd43a9abc1f2f18e08d391\c11c8acdc82b4f71b2eb4c4ab3a15496\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite
ERROR: Raster map or group <output18d498a864564b1ebfcce3671ea42038> not found
C:\PROGRA~1\QGIS3~1.4\bin>exit
Execution of <C:\Users\seansuzedielois\AppData\Local\Temp\processing_5a9b7eb19ebd43a9abc1f2f18e08d391\grassdata\grass_batch_job.cmd> finished.
Cleaning up temporary files...
Press any key to continue . . .
Execution completed in 1.60 seconds
Results:
{'output': <QgsProcessingOutputLayerDefinition {'sink':C:/Users/seansuzedielois/AppData/Local/Temp/processing_5a9b7eb19ebd43a9abc1f2f18e08d391/c11c8acdc82b4f71b2eb4c4ab3a15496/output.tif, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
The following layers were not correctly generated.

C:/Users/seansuzedielois/AppData/Local/Temp/processing_5a9b7eb19ebd43a9abc1f2f18e08d391/c11c8acdc82b4f71b2eb4c4ab3a15496/output.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts:
i am having some trouble and i could use some help. I downloaded a Geotiff file on WGS84 for the area that i am working for and i would like to clip it so it fits the outline of my map (red polygon).

All the layers of my map are set on Greek Grid. I use the Project tool from arc toolbox (Projections and Transformations-Raster) and i create a new shapefile with GreekGrid as spacial reference. Then i clip my file with Management Tools -> Raster -> Raster Processing -> clip.

This is the result i am getting. The red box is the outline of my map and it doesn't align with my file.

Am i doing something wrong? Someone suggested to use the extract by mask by mask tool because but i am having odd results with this tool. I post screen shot of the result below below. The first on is with using the clip tool, and the second with using the extract by mask tool.


texts: Set projection/CRS in map and shapefile to produce buffer

texts:
I simply want to create a buffer of a specific size, around a single point. When I create a buffer it uses degrees/ or radians, or whatever, and the buffer is created that covers the whole screen.

I am working on a South African map, so I have tried setting the project CRS to WGS84/UTM30S.

I then load my shapefile, and set the shapefile CRS to WGS84/UTM30S. Still the buffering tool creates this massive buffer.

I just want to be able to create a buffer, of known radius in meters.

I am using QGIS 1.8.0 Lisboa

Anyone with ideas?


texts: Modifying "Multiple Ring Buffer" Magnitude based on value of input point

texts:
I'm trying to visualize the area of influence around an earthquake epicenter based on Magnitude.
Is there a way to modify the multiple ring buffer function so that it's not just a given distance away from the epicenter?


texts: GeoTIFF gets pixelated when rendering on GeoServer?

texts:
I am trying to render an 8 bit GeoTIFF (4 Band) on GeoServer and I am having troubles rendering the image on an WMS service.

When I zoom out the image, it gets pixelated:

Below, it is the same GeoTiff in QGIS with same scale (1:100000)

The image rendered by GeoServer is extremely pixelated. I tried tweaking the configuration but I could not reach the same level of quality as it is shown in QGIS.

The version of GeoServer we are using is the 2.2.4. Could someone give me a hint about what may be causing such pixelation on GeoServer?


texts: Customizing map layers in Esri Javascript

texts:
I'm working on creating an outdoors-app-thingy where I'm required to divide up parts of 3 different rivers where the app is focused around into different "zones".

After a lot of calling around and searching, I came across ArcGIS/ESRI which has everything that I need and it's a lot more detailed than anything else out there, and have been playing around with it for a few days now.

I was wondering if it's possible to add something similar to click events to the Areas that I'm adding to the map?

The thought behind it is that the page you navigate to the map from transmits the data that sets the map center and zoom, depending on which river you selected on the previous page you have 1 to 3 "areas" that are drawn on the map.

When clicking one of these "Areas", jQuery would simply just direct you to the final view where you'd get more streamlined information.

Would something like this even be possible?


texts: Using QGIS layers in ArcGIS SDK for Android

texts: So I use QGIS for work and I am trying to create app for field work. I'm looking into ArcGIS SDK for android, but it looks like you will need to import your layer into ArcGIS server. Is there any NON-ArcGIS service I can use to import my layer onto my app, still using the sdk.


texts: Reclassifying multiple rasters in ArcGIS Desktop?

texts:
I have 163 rasters, each one representing the least-cost path from a source to multiple destination locations. I would like to reclassify each raster in order for the cells representing the path to have value 1. I want to do that since I am planning to sum all the rasters up and to use a color scale that would make evident (i.e., more intense color) where paths overlap.

Now, the problem is that I do not know how to reclassify the rasters in a single operation. I am wondering:

(A) can the "iterate rasters" function in the ModelBuilder facility be used for my purpose (and if it can, how?);

(B) how can the Reclassify tools in Arcgis 10.1 be put to work in conjunction with "iterate raster"?

Furthermore, I anticipate that an additional issue would be that the rasters have not the same values. In some of them, the values span from 1 to 32, in other ones from 2 to 28.


texts:
It seems I have a similar problem to GeoTIFF world file ignored by QGIS?

QGIS doesn't import my ground-layer tiff-files correctly, because I have the geo data in an extra tfw-file.

Now I want to delete all gdal-metadata in a bunch of tiff files.

Is there a Linux command line tool for this?


texts:
I have raw .dat file which i need to convert to a raster.
I am currently using numpy to read the file
data = np.fromfile('gsmap_gauge_now.20200805.0400.dat',dtype = np.float32)

The coverage of this raster is from 60 N - 60 S with grid of [3600 x 1200]
Can anyone please help me out on how to achieve this using GDAL ? or any other linux supported tool ?
Thanks a lot


texts: Google Earth Engine; Upload TFRecord via https://code.earthengine.google.com/ or command line fails

texts:
The Files: Json TFrecord

The TFrecord is the predict result from a tensorflow model. Its a [1,1] Numpy array per patch and has been transformed to TFRecord by(like in this example) :

`for predictionPatch in predictions:
    print('Writing patch ' + str(patches) + '...')
    predictionPatch = predictionPatch[
        x_buffer:x_buffer+KERNEL_SIZE, y_buffer:y_buffer+KERNEL_SIZE]

predictionPatch = predictionPatch[
    x_buffer:1, y_buffer:1]
# Create an example.
#example = tf.data.Dataset.from_tensor_slices(predictionPatch).SerializeToString()

example = tf.train.Example(
  features=tf.train.Features(
    feature={
      'OC': tf.train.Feature(
          float_list=tf.train.FloatList(
              value=predictionPatch.flatten()))
    }
  )
).SerializeToString()

print("example",example)
print("type(example)",type(example))
#print("type(tf.io.parse_single_example(example))",type(tf.io.parse_single_example(example)))
# Write the example.

writer.write(example)
patches += 1

writer.close()`

uploading the file via browser mask https://code.earthengine.google.com/ works

but the ingestion results in :

changing the shape of the base date to 1 and :

predictionPatch = predictionPatch[
    x_buffer:1, y_buffer:1]

to :

predictionPatch = predictionPatch[
    0]

results in a different error in the browser mask:

Manifests for TfRecord ingestion must have exactly one tileset with exactly one source.

Uploading via earth engine on Linux:

earthengine upload image --asset_id=users/AVHomefolder/XCeption_Carbon_KS1 gs://masterarbeit_bucket/test/XCeption_Carbon_KS1.TFRecord gs://masterarbeit_bucket/test/XCeption_Carbon_KS1.json

(the flag --no-use_cloud_api will start but will not upload)

throws some Dependencie Errors, but starts the Job.
this error is retrieved from the browser mask:

reloading the TFRecord and using .take(1)[0] results in :

<TakeDataset shapes: {OC: (1, 1)}, types: {OC: tf.float32}>
{'OC': <tf.Tensor: id=31, shape=(1, 1), dtype=float32, numpy=array([[-0.24551529]], dtype=float32)>}

using the same on the equivalent file from the example results in similar Tensor :

{'OC': <tf.Tensor: id=31, shape=(1, 1), dtype=float32, numpy=array([[0.]], dtype=float32)>}


texts: Converting zoom level rule-based styles to mbtiles?

texts:
In QGIS I have made a zoom level rule-based style (.qml) with OpenStreetMap data stored locally for a specific area.

Is it possible that I convert these styles (or the project) as mbtiles so that I can host it as a basemap for my Leaflet project?


texts: Making Layers change color based on attribute expression in GeoServer?

texts:
At the moment I'm working with a dataset that has a count value as the only value in the database.

However, if I add &VIEWPARAMS:Sex=Male to the given url from GeoServer, the results change based on this query.

I want to change the layer of GeoServer based on what is given in this query.

For instance, if Sex=Male I want the color of the layer to be blue and if Sex=Female I want the color to be pink (in order to give a beter representation and eventually split up the rasters).

Does anyone have a solution for this, because I have been looking for days and I have not found any possible solutions?

-EDIT-

Added both pictures the picture in pink is how the dataset for the female set should look like and the blue one is the male set. What I now want to realise is that if I change the url by adding: 'VIEWPARAMS=Sex:male' I want the color of the map to turn blue and if I do 'VIEWPARAMS=Sex:female' I want the color of the map to turn pink. Without adding another kind of command like '&layers=Male_Layer' or '&layers=Female_Layer'. Is this possible in geoserver?

Incase the pictures are not working here are the url's http://imgur.com/VCrkoss and http://imgur.com/bmFNFCB


texts: Store in dictionary/list strings with special characters from Attribute Table

texts:
I need to store into a dictionary data from a shapefile that contains in several fields unrecognizable characters by ASCII, like in this one: "Städtebaulicher Entwicklungsbereich: Sanierungsgebiete SG 1 u. SG 2, Mkt. Floß". In the end, I have to create a report in Excel with all these attributes.

I get the data ( stringVal ) by using SearchCursor. This is a simplified code of how I store the data in the dictionary:

cursor=arcpy.da.SearchCursor(filename_path,["Category"])
for row in cursor:
        stringVal=row[i]
        if stringVal not in d:
            d[stringVal]=1
        else:
            d[stringVal]=d[stringVal]+1

So my dictionary should look like this:

d={"Städtebaulicher Entwicklungsbereich: Sanierungsgebiete SG 1 u. SG 2, Mkt. Floß":1,"Städtebaulicher Entwicklungsbereich: Sanierunggebiet Südstadt, St. Grafenwöhr":4}

At this moment, it looks like this, when printing it:

d={u'Bayerisches Staatsministerium der Finanzen, f\xfcr Landesentwicklung und Heimat': 11}

I know it has to do with unicodes and UTF-8, but from what I read on the internet, I still couldn't solve it.


texts: Merging one DEM into another to fill in missing data using QGIS?

texts:
I have a LiDAR dataset with some missing data right in the middle of my map. I would like to use a 5m DEM to fill in the missing areas.

Is there a way in QGIS to merge the two but only use the 5m data in the missing location?


texts: Show imported waypoints in table of contents (QGIS)

texts:
I collect data with a garmin handheld GPS, and then import said data into GIS software. For a specific project, I need to manipulate the buffers for individual points that Ive collected.

Is there a way to showcase each point in the table of contents? For now, they all show up as one layer... but I would like to express different symbologies for each one.


texts:
I have 397 shapefiles that I need to merge into one shapefile, and I need to add an attribute for the name of each of the shapefiles. I have no experience in coding/ Python.

Is this possible to do without any coding in ArcMap 10.3?


texts: DataReviewer Geometry check on String and Int

texts:
I just started using Data Reviewer and went to compare two feature classes. A polygon of Parcels and a feature class of Points of addresses. I want to compare the address to ensure the address number is the same where they are on top of each other. The problem is that the Address field in the Parcels Feature Class is a string whereas the Address Feature class is an integer. I need to be able to convert the Addresses field in the points feature class to String just for this check.

Let me know if you need anymore information.

Edit 1: For the parcels feature class i have the address column with values such as 1164, 1160, Null, 1-19/1-26/3-1 , -/-, -/0/- (those are all the types total of 1062 values)

Then there is Address_number in the points feature class which is a long. filled with: Null and 1,3,542,6785,1234 etc. (all the types.).


texts: Calculating area of vegetation in raster file?

texts: How do I calculate area of several vegetation types in a raster (.img file type) image. They are color coded. Am using QGIS (1.8.0) and ERDAS 9.2.


texts: Warning message using writeRaster() - not updated for proj >= 6

texts:
Using the code below to write a raster to disk:
library('raster')

r <- raster(nrow = 5856, ncol = 9141)
r[] <- 1:ncell(r)

writeRaster(r, tempfile(fileext = ".tif"), overwrite = TRUE)

always produces a warning message:
In .gd_SetProject(object, ...) : NOT UPDATED FOR PROJ >= 6

What causes the message and how can I avoid it?


texts: How to Select all polygons or lines that lie at least partially within a circle with a radius?

texts:
I am storing OS Open Local data, as downloaded here: https://www.ordnancesurvey.co.uk/opendatadownload/products.html (top result).

In particular my geometry column takes one of two possible shapes:

SELECT ST_AsGeoJSON(geometry) FROM roads LIMIT 1;

Resulting in:

{"type":"LineString","coordinates":[[464223.46,1213396.46],[464257.63,1213398.4],[464294.12,1213438.78],[464330.93,1213444.89]]}

And

 SELECT ST_AsGeoJSON(geometry) FROM buildings LIMIT 1;

Resulting in:

   {"type":"Polygon","coordinates":[[[454003.53,1203047.09],[453984.86,1203083.28],[453997.27,1203089.71],[454015.96,1203053.49],[454003.53,1203047.09]]]}

Given a lat, a long, and a radius in meters, how would I SELECT all roads or buildings which are either contained by, or intersect the selected circle?


texts: Affecting graphics in MXD by geographic transformation using ArcGIS Desktop?

texts:
Our GIS users have been using graphics to mark things in their map documents. By graphics, I mean circles, rectangles, lines, text, etc. that can be drawn into data view by using the Draw toolbar. We have a spatial reference (sr) change of all our geodata upcoming. Using Python, I want to crawl through hundreds of MXDs and replace all layer's data sources to point to the feature classes in the new sr. Also the script will set the sr of the data frames to the new sr.

But what about the graphics?

After the sr change, the graphics must still lie correctly relative to the features in the new sr. I didn't find a way to directly access and manipulate graphics with arcpy. As a workaround, I tried creating a custom transformation "sr old to sr new" and added it to the dataframe, assuming that the graphics have been drawn in the old sr and should therefore be transformed to the new sr. Unfortunately, this had no effect on the graphics. After many test runs of the script, the graphics lie displaced relative to the features in the new sr.

How can I transform the graphics to the new sr, so that they will maintain their correct position relative to the features?

Is there a way to directly access and manipulate graphics with ArcPy?

I am using ArcGIS Desktop 10.2.1.


texts: How to text in selected polygons using QGIS in geojson file?

texts:
I have upload geojson file in QGIS. there total 8 polygons (8 taluka's co-ordinates) now, my question is I want to just highlight 4 polygons( 4 taluka's) and display in selected 4 polygons ( "number of villages " &" area for crop" ) for selected 4 crops . and display names in the each polygons also. is it possible ?

I want to display  in "Kannad" :  "Number of villages :  12, "Area for Crop" : 20 Hectare "  only


texts:
I am trying to create a shared server for community, and now stock=k with users privilege and roles, I get an error on QGIS when login with my user account.

Steps I did :

I revoked all privileges from public schema.

REVOKE ALL ON DATABASE sudanmap FROM PUBLIC;
REVOKE ALL ON SCHEMA public FROM PUBLIC;
REVOKE ALL ON ALL TABLES IN SCHEMA public FROM PUBLIC;
REVOKE EXECUTE ON ALL FUNCTIONS IN SCHEMA public FROM PUBLIC;
I created role named userread only to view as demo.

CREATE ROLE userread;
Granted Connection to database, usage to public schema and Grant Select .

GRANT CONNECT ON DATABASE sudanmap TO userread;
Grant USAGE ON SCHEMA public TO userread;
GRANT SELECT ON ALL TABLES IN SCHEMA Public TO userread;
Created user with password, and grant userread to it.

CREATE USER open PASSWORD '123456789';
GRANT userread TO open;
Login to QGIS is successful. I can see my schema but I get this (failed to get layer).

What I try I grant all to my user still the same, only why he can see the layers by making him  superuser.

What am doing wrong ?

Edit: All Sequences grant to user, but still the same problem


texts: Getting multiple value pairs as Python Script tool parameter?

texts:
I've tried several different data types but can't figure out how to allow a user to create a list of elevation-temperature pairs. For example, a spreadsheet might look like this:

Is there a parameter data type that would allow the user to specify these values as pairs?


texts: Export from a GIS database and automatically populate an external document

texts:
I'm using a standard database based on ESRI polygon shapefile with 30 fields. The database is diverse with records from text values to long integer and double numbers. The records from the database needs to be linked with an external document/smart pdf file or Word file.  The Word file/smart pdf is formatted in a manner that only some parts of the document needs to be completed.
For example, I need to complete/autocomplete a sentence or a table with records from ESRI polygon shapefile, such as : Name ..... , Streets ......, city ........

Is there a plug-in, tool or algorithm scheme in order to link records from a shapefile (in software such as ESRI ArcGIS or Quantum GIS) to particularly paragraphs from an external document?


texts: Showing lat long instead of numbered index on spatial data map

texts:
I am trying to project meteorological data from a netcdf4 file with its latitude and longitude information and I am not getting how to show lat-long information on the x and y-axis, in the image instead of index values.

Here is the code I have written:
import gdal
import matplotlib.pyplot as plt
# Path of netCDF file
netcdf_fname = "F:\TS\d02_2021-03-31_000000.nc"
# Specify the layer name to read
layer_name = "mcape"
# Open netcdf file.nc with gdal
ds = gdal.Open("NETCDF:{0}:{1}".format(netcdf_fname, layer_name))
# Read full data from netcdf
data = ds.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)
data[data < 0] = 0

A = data[4,:,:]
D = plt.imshow(A)
plt.colorbar(D)
plt.title('mcape',fontweight ="bold")
plt.show()


texts: Arcpy unzipping script: ERROR 000714 Error in executing: cmd.exe /C

texts:
I wrote script for unzip archive with GDB.
It works locally, i.e from ArcGis Desktop.

But when I publish scrypt as geoprocessing service, it doesn't work - ERROR 000714: Error in script Zip22. Error in executing: cmd.exe /C

Script:

import arcpy
import os
import sys
import traceback
import zipfile
import re

class LicenseError(Exception):
    pass

def unzip(zipFileName):
    try:
        path = os.path.dirname(zipFileName)
        zip = zipfile.ZipFile(zipFileName, 'r')
        gdb_name = None
        for name in zip.namelist():
            if '.gdb' in name.lower():
                gdb_name = name.replace('/', os.sep)
                break
        if not gdb_name:
            raise Exception('zip file has no gdb')
        zip.extractall(path)
        zip.close()
        return path + os.sep + gdb_name
    except RuntimeError:
        zip.close()
        arcpy.AddWarning(get_ID_message(86133))

def get_ID_message(ID):
    return re.sub("%1|%2", "%s", arcpy.GetIDMessage(ID))

if __name__ == '__main__':
    try:
        zipFile = arcpy.GetParameterAsText(0)
        gdb_name = unzip(zipFile)
        arcpy.SetParameterAsText(1, gdb_name)

    except:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        pymsg = "ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n        " + str(sys.exc_type) + ": " + str(sys.exc_value) + "\n"
        arcpy.AddError(pymsg)

Similar script with archiving GDB to zip works fine both locally and at server.

What's wrong in this script?


texts: Check if point is inside unclosed linestring using POSTGIS ST Functions

texts:
I have a table with spatial data. This table contains a single record which is a unclosed linestring.

I calculated the Centroid of this geometry using query

SELECT ST_AsText(ST_Centroid("geom")) FROM "public"."ha_ground"

I got the output "POINT(73.8180699704927 18.6281375643436)"

My aim is to check if this point is inside the unclosed linestring

I used the following query

SELECT ST_Contains("geom",ST_GeomFromText('POINT(73.8180699704927 18.6281375643436)', 4326)) FROM "public"."ha_ground"

As We know Centroid of this geometry is inside the unclosed linestring.
But, This query returns false instead of true.

What is the possible reason for this.

And what is the solution to check if point is inside unclosed linestring


texts: How to join a buffer (point) to a line in QGIS 3.0 modeler?

texts:
Is there a way to join an output buffer of a point to a line in QGIS 3.0 model processor? This is being added to a larger model so a standalone script wouldn't be practical.

I am not sure if this is an error that needs to be reported or if I am doing something incorrectly.

I am wanting to buffer a point and then join it to a line in QGIS 3.0 processing modeler. However I cannot set the output of the buffer as the input or join layer for the "Join Attributes by location" algorithm.

I have attempted the following things:

(1) Doing each step individually outside of processing modeler (works)

(2) Entering the input myself

(3) trying the same with a single sided buffer (works)

(4) selecting the buffers by location and then joining (I am able to select, but then not join).

(5) joining the points to transmedia before buffer (works)


texts:
I am pretty new to GIS/Mapping apps. I have obtained some parcel data/tax data and imported it into SQL using Shape2SQL.

The values that I have in the DB are based of the following projection definition.

PROJCS["NAD_1983_StatePlane_Florida_West_FIPS_0902_Feet",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",656166.6666666665],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-82.0],
PARAMETER["Scale_Factor",
0.9999411764705882],
PARAMETER["Latitude_Of_Origin",24.33333333333333],UNIT["Foot_US",0.3048006096012192]]

The sample data has values such as:

POLYGON ((713262.40299399942 725948.04824899882, 713242.40249400213 725948.04974899441, 713242.36924400181 725908.74974899739, 713291.95949400216 725908.74599899352, 713291.99174400046 725933.04624899477, 713262.40149400011 725933.04824899882, 713262.40299399942 725948.04824899882))

I get these values by using the geom.STAsText() function.

I am trying to use Bing Maps, which passes in a polygon of Lat/Long Coordinates that represents the bounds of the map. I want to show parcels that are within the bounds of that view (POLYGON ((-82.005076780915289 26.604321062391364, -82.005076780915289 26.586956646764229, -81.964822188019781 26.586956646764229, -81.964822188019781 26.604321062391364, -82.005076780915289 26.604321062391364)))

I am really confused about data types and I can't seem to pull and results from the DB. I am thinking the False Easting is also something I need to take into account.


texts: ArcGIS color bar from pre-specified RGB values

texts:
I have a 256 item vector of RGB triplets derived from viscm (a tool for creating perceptual color bars in a CIE color space). I am not able to replicate it correctly using the algorithmic color bar facilities in ArcGIS, even when using the CIE LCh and LAB options (because the path taken through the color space is an arc rather than a line, as seems to be assumed in ArcGIS).

Is there a way to specify a color bar by giving it my 256 values directly? If I wanted a multi-part color bar, is there a way to specify that in text format for so that I can interpolate between each of my 256 colors, without having to do it through the GUI?

As a first step, this is designed for raster symbology, but should also be applicable to vectors. Note also that I don't want a discrete colour map, but one that interpolates between the values I specify, ideally in one of the CIE spaces.


texts: File iterator in ModelBuilder does not work correctly when passing variable input folder?

texts:
In ArcGIS ModelBuilder, I am trying to iterate over some ascii files (my_ascii.asc) in a folder (c:\my_folder) and converting them from ascii to raster. This works fine if I do it as shown below.

however I want to pass some parts of the input folder as a model parameter.

So let's say my input folders are c:\my_folder_2010 and c:\my_folder_2030, i want the year to be a model parameter. If I make the year variable and pass it into the iterator like c:\my_folder_%year% as shown below, it does something weird. Iterator finds the correct files, outputs the correct 'name' ('my_ascii'), but the output ascii remains 'file.asc' instead of 'my_ascii.asc'. Once it enters the 'ascii to raster', it says that it cannot find the 'file.asc', which is correct, because it is called 'my_ascii.asc'.

Why does iterator only works correct for the 'name' part, and not for the 'file.asc'? Any solutions?

I also tried parsing a path into the first data element that goes into the iterator, but that also did not work.


texts: How to convert fake donut holes in polygons to real ones?

texts:
I have this problem where I have to import shapefiles into our system which must be geometrically valid, but the software that everyone uses to produce them seems to output files where every polygon that should contain a hole is instead made up from a single outer ring which goes around the boundary and then goes inward along some arbitrary path, traces around what should be the inner ring, and then doubles back on itself back to the outside.

What I seek is a way to resolve these by somehow dissolving the vertices along the arbitrary path in and back out of the hole, to make a seperate inner and outer ring which will form a valid polygon. I have tried the trial version of ArcGIS and it seems to do this automagically when you snap the vertices on the double-back part together, but my business doesn't have a copy of ArcGIS and I might have a hard time convincing my boss to buy it, so I was hoping QGIS would allow me to do such a thing?

Screenshot:

I couldn't figure out how to get any WKT out of this sorry, but from playing around with the vertices I'm quite sure that it's one single ring that doubles back on itself for a bit.


texts: Why does shp2pgsql generate a MultiPolygon from a Polygon shapefile?

texts: I created a shapefile with a single Polygon in it using QGIS. When I test it with ogrinfo it says it is of type "Polygon". When I convert to GeoJSON using ogr2ogr it also comes out as Polygon. But when I use shp2pgsql it comes out as a MultiPolygon. What's the reason for this?


texts: I am using QGIS 2.14.1 in Ubuntu 15.04 and when i click Raster -> HeatMap -> HeatMap. My QGIS crashes and closes automatically even there is a point layer loaded in the canvas(HeatMap Plugin expects a point layer). Any solutions for the same?


texts: PostGIS clipping raster with the same color band as the original raster

texts:
I am trying to clip a raster so that the clipped raster would have the same minimum and maximum values for the color band as the original raster.

For example if i clip the raster with this command:

create table clipped AS SELECT ST_Clip(rast,ST_MakePolygon(ST_GeomFromText ('LINESTRING(424920 7370964, 424920 7371476,425432 7371476,425432  7370964,424920 7370964) ')) ) FROM original;

I get the following result:

So the clipped raster has band minimum value as 0 and maximum value as 20. Instead of that, i want it to be exactly same as the 'original' raster: from 0 to 33. In other words, I want the color of the clipped raster to be same as it is in that area in the 'original' raster.

I am new for PostGIS, but i have tried to find an answer for a while to this question but i haven't found anything. I have tried the functions like ST_Reclass but with no success. What is the proper way to do this with PostGIS?


texts:
I want to create hotlinks to open excel and image files but am having issues.

I fixed the issue with opening excel files by adding ""

cmd /c start excel.exe /r  "[% "Excel" %]"

However I can't get images to open??? rather than opening a specific program can qgis just open a file with the default system viewer?

cmd /c start rundll32 "C:\Program Files\Windows Photo Viewer\PhotoViewer.dll" "[% "Location" %]"

I tried a few option with rundll32.exe C:\WINDOWS\System32\shimgvw.dll,ImageView_Fullscreen %1 but they don't work.

Also tried  as per Answer from Nathan but I get "Unable to run command"

=== Original post follows.

I have looked at Adding a hotlink to an attribute field in QGIS and http://osgeo-org.1560.x6.nabble.com/Actions-HowTo-td4534256.html . The code from here launches excel but can't find the file...

I have the following in the properties-actions

when I try to open it I get

It looks like the issue is with the spaces in the directory name...

P:\2013\361_xxx_Production\Input\Supply1\28-05-xxx_Application_Packages\Roads photos and form 2s\Roads Package 4\Form 2 Output xxx Roads Package 4.xls

it seems to be getting truncated at the first space in the directory. It is bad datastructure but I can't change this. Is there any way to fix this issue?


texts: Making colours from table with colour percentages/RGB-codes

texts:
Currently, we have a big dataset imported in our QGIS profile. Normally we would have imported this by WMS but because this particular WMS/WFS consists of more than 20 parts and the data is technically complete and will never change. We could luckily download the vector data directly. It is called the "Tertiair geologische kaart".
Which can be found here: http://www.geopunt.be/catalogus/datasetfolder/06d1a1d7-d46c-43d6-96aa-25225b57dd8b.
But here comes our problem, for this dataset no SLD or similar kind of file was constructed to make up the colour of the legend. We know how the legend is constructed on a specific field, **so that is not the problem is of the colour, which has been set as a standard for this dataset.
Normally you get this with the WMS/WFS but that is not the option here. The legend is also supposed to be "Listhographic".
How do we set the symbology colour based on the colour percentage/RGB table?
There is a table that does have the colour percentages and RGB-codes, we managed to link this table based on the "code column". So, in theory, we should be able to just set symbology based on those columns and values, while still maintaining the legend as it is but the problem is we do not how to do that.
Here is a sample of the table:

Here is a sample of the legend:


texts: How to "join" or "merge" two polygons in QGIS and fix "overlapping" points (green X)

texts:
I'm trying to digitize polygons from a map and I now decided to "join" two polygons that were separated originally. Here is a picture showing a polygon with the line which separates the two polygons (same feature). I want to "join" "merge" or "dissolve" the polygons but I wasn't able at all to do it. How can I perform a dissolve on this feature so that I get the proper joined polygon? How can I do this without creating a new layer and "overwriting" the layer in which the polygon is found?

While trying, I ran into a problem saying my geometry is not working. But I'm not able to get rid of the green X where three points meet creating this problem. Below is the faulty point. How can I correct the three overlapping points so that the geometry is fine?


texts: How to display point features w.r.t field values of a polygon file in arcobjects?

texts:
I have a point & polygon shapefile. Point file has a field namely, Type: atm,hospital,school etc... Polygon file has a field Zone: Zn1, Zn2, Zn3. There will be two combo boxes which will populate this two fields respectively

What i like to do is.......If i select school & Zn2 in the comboboxes and click the "Show" button...... only schools which fall under the Zn2 polygon should be displayed. W.r.t arcobjects, i m just a beginner. Trying this for the past two days. Can anyone help me regarding this ? Thanks in advance...

(ArcGIS Version 9.3 )


texts: Asking for the precise detail of gdal_translate's output .xyz multiband file format (when converting from GRIB)

texts:
gdal_translate is able to convert a multidimensional GRIB data file (Copernicus ERA5 reanalysis) with the following command-line command (gdal_translate.exe from 11/2016):
gdal_translate source.grib target1.xyz

gdal_translate will in this case create two files:

target1.xyz.aux.xml
target1.xyz

target1.xyz.aux.xml is rather clear, it contains readable data like this:
<PAMDataset>
  <PAMRasterBand band="1">
    <Metadata>
      <MDI key="GRIB_COMMENT">10 metre u wind component [m/s]</MDI>
      <MDI key="GRIB_ELEMENT">10U</MDI>
      <MDI key="GRIB_FORECAST_SECONDS">0 sec</MDI>
      <MDI key="GRIB_REF_TIME">310262400 sec UTC</MDI>
      <MDI key="GRIB_SHORT_NAME">0-SFC</MDI>
      <MDI key="GRIB_UNIT">[m/s]</MDI>
      <MDI key="GRIB_VALID_TIME">310262400 sec UTC</MDI>
    </Metadata>
  </PAMRasterBand>
  <PAMRasterBand band="2">
    etc.

target1.xyz is a mystery for me. It's definitely not readable ASCII, but binary. It's not possible to paste the first 100 bytes into this SE edit window. Reading from file start as 64-bit doubles produces crazy values (huge exponents), so it definitely has a structure and most likely a header.
No solution:
Once I have the xml file, I can parse it's headers and get the corresponding band numbers. Then I can call gdal_translate again, this time with a band number (note "target2" file name). Example:
gdal_translate source.grib target2.xyz -b 12345

In this case, I do get readable ASCII, example:
22    60     4.953857421875
22.25 60     5.292724609375
22.5  60     5.439208984375
22.75 60     5.296630859375
23    60     5.023193359375
23.25 60     4.396240234375
23.5  60     3.891357421875
23.75 60     4.203857421875
24    60     4.513427734375
24.25 60     4.732177734375
24.5  60     4.951904296875
etc.

However, each conversion takes one second and there can well be 30000 bands to extract => 8 hours! I would very much like to read the target1.xyz file using own code, but it seems impossible to find any information about it's format.
Tried in the first conversion ("target1", with no band argument):

Using other target1 file name extension, in order to get another output format - gdal_translate ignores.
Specifying the output format with the -of [format] parameter, gdal_translate ignores.
Reanaming target1.xyz => target1.tiff and loading into a program with capabilities - no go (says invalid file format). Seemed like a waste of time to try more renaming + loading, blind shooting.

Tried in the second conversion ("target2", with band argument):

Specifying multiple bands simultaneously (as the specification says it is possible, see
https://gdal.org/programs/gdal_translate.html) - gdal_translate ignores all but the first band (this is also verified elsewhere in SE); in addition the DOS command line has a length limit. (I could have imagined that gdal_translate had created multiple output file name variations in this case). Example:
gdal_translate source.grib target2.xyz -b 12345 12346 12347 etc.

Question: What is the file format and structure of the target1.xyz file that
gdal_translate source.grib target1.xyz

creates?