texts: Use of Leaflet plugins in 'Leaflet for R' to connect with ArcGIS online

texts:
I was wondering if it is possible to use the plugins for the normal (non-R) version of Leaflet, and if, how I can do that.
reason I need to do that is, that some of the data I need to show are hosted at ArcGIS Online.

Currently available plugins via leaflet.com are 'leaflet-providers' and 'Esri Leaflet'


texts: Using QGIS GUI predefined colorramps in Python console / PyQGIS?

texts:
When using the QGIS Gui, there are plenty of predefined colorramps to choose from when styling a layer with e.g. categorized symbology.

When using the Python console for PyQGIS scripting, this page of the API docs describes how to create your own colorramp that you can then apply to a renderer object, for example:

# create a new colorramp
colorRamp = QgsVectorGradientColorRampV2.create({'color1' : '0,255,0,255',
                                                     'color2' : '255,0,0,255',
                                                     'stops' : '0.5;255,255,0,255'})

# set up an empty categorized renderer and assign the colorramp
renderer = QgsCategorizedSymbolRendererV2(field, [])
renderer.setSourceColorRamp(colorRamp)

Is there a way to access the predefined colorramps that are accessible from the GUI and assign them to a renderer instead of creating my own?


texts: Getting data from Mapire.eu (website) to QGIS

texts:
Is it possible to get data from this web site (mapire.eu) to QGIS?

I want to make a historical map of rivers in 18th and 19th century compared to today but cant find a way to import data.

I am new to QGIS.


texts: Choosing a GIS which imports a list of stage addresses, calculates the route, allows to modify and export it

texts:
I'm new to GISs. I'm searching for a GIS to meet following requirements.

I have a list of street addresses [A,B,C,D,E,...], with additional information associated (basically which orders to deliver at each address and to which customer), available in a DB table (DBMS: DB2 for i Series) or .xls file.

I need to show these addresses on a map and calculate the route among them. Not the global best route, but the route that touches each point maintaining the order given for them. Then the user should modify this order if necessary (e.g.  [A,E,C,B,D,...]) and the route should be recalculated.

As output I need to write the points in the new order on a DB table (also with customer info etc. for each record). I also need to print a PDF with the map and route instructions for the truck driver.

I considered MapPoint, but counld't succedd with automatic route calculation for the imported addresses (to calculate a route I have to manually add each point). Google Maps API seems promising, but I don't know how to do many of the operations above.

Any hint on which GIS could fit? Both web-based or desktop applications are welcome. Many thanks!


texts: How to convert a shapefile to a grid file

texts: I was wondering whether anyone might be able to tell me in layman's terms how to convert a polygon layer to a Grid file in ArcGIS 10.1 or ArcGIS 10.2?


texts: Migration from postgresql 8.4 to postgresql 9.0, postgis 1.5 involved

texts:
Im migrating from 8.4 to 9.0, in pg8.4 i have postgis 1.5.1 and in 9.0 im just installed 1.5.2.

I decided to dump the postgis databases before the pg_upgrade command, and delete them, planning to create my new postgis_template and then create the databases:

pg_dump -C -h host -U postgres -f databasegis.dump databasegis
psql -U postgres < databasegis.dump

But I have errors :

ERROR:  could not access file "$libdir/liblwgeom": No such file or directory
ERROR:  function public.st_text(boolean) does not exist

What can i do?

Thanks a lot.


texts: Numbering rows and columns in polygon grid (fishnet) using ArcGIS Desktop?

texts:
I seem to be a bit stuck on what I thought would be quite straightforward.

I've created a fishnet of polygons, all 25 x 25 metres, 100 polygons by 100 polygons.

My main task is to create a relative co-ordinate system based on the distance from a vertex somewhere near the middle of the grid (but NOT the middle), which will be the origin 0,0.

Basically, the polygon square - where the origin sits at the very bottom left of it - will be X=0 and Y=0. The next square to its right will be X=25 Y=0 and so on and so on to the extents in all directions. The increasing blocks of 25 just represent the actual distance in metres from the origin.

I thought I could maybe number the columns and rows of the polygon grid then work from there, assigning an X and Y value based on column and row number but I can't even seem to do that.

Any ideas?


texts: Strange result with ST_Distance, geography type and equator

texts:
I am getting strange results with distance queries in PostGIS using geography types and ST_Distance. According to the docs, ST_Distance should return the distance in metres between two geography objects. However this query, between a polygon and a point that is two degrees south of the polygon, returns a distance of 0:

SELECT ST_Distance(
    ST_GeographyFromText('POINT(0 -82)'),
    ST_GeographyFromText('POLYGON((-90 -80, 90 -80, 90 10, -90 10, -90 -80))')
)

I've verified that the lat/lon order for the WKT strings is as expected (lon/lat order; x/y). The polygon is a valid geometry. Both ST_Touches and ST_Intersects return false between this point and polygon geographies, which suggests to me that the distance must be greater than 0. But the result is stubbornly 0...

That is until you change the polygon to not cross the equator (no, I don't mean the antimeridian). If you change the northern edge of the polygon to be south of the equator like below, the result of the ST_Distance query is 893,493.4 metres, when it shouldn't have changed at all.

SELECT ST_Distance(
    ST_GeographyFromText('POINT(0 -82)'),
    ST_GeographyFromText('POLYGON((-90 -80, 90 -80, 90 -10, -90 -10, -90 -80))')
)

In either case, if I cast to geometry, I get 2 (degrees).

This really stumps me. I can consistently reproduce this result in two different environments:

"POSTGIS="2.2.1 r14555" GEOS="3.5.1-CAPI-1.9.1 r4246" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 2.1.2, released 2016/10/24" LIBXML="2.9.3" LIBJSON="0.11.99" RASTER"
POSTGIS="2.3.0 r15146" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" TOPOLOGY RASTER


texts: How to force ol3 map to immediately reload tiles from geoserver after a change?

texts:
I am having a trouble with the cache of the browser, when I make changes, they don't appear until the next animation (ie zoom change ). What I want is an immediate reload of images after changes right now am using this method:

layer.getSource().updateParams({
  'ol3_salt': Math.random()
});
layer.getSource().setTileLoadFunction(layerZdcIntWMS.getSource().getTileLoadFunction());
layer.getSource().changed();
map.updateSize();

but it doesn't reload immediately
what should i do ?


texts: Migrating ArcView GIS *.apr to ArcMap application of ArcGIS Desktop 9.3?

texts: I have a lot of arcview .apr file. I need to open & edit in ArcMap 9.3, but the original project file is saved in .apr. I tried to use "import from arcview project". The project file can be open finally, but some data source was broken. I cannot trace the original data file, as the layer name was assigned by anothers. How can I trace the data file and re-link it in Arcmap 9.3?


texts: Print Composer exporting different image with Thunderforest base map

texts:
I prepared a map in the Print Composer of QGIS. The scale of the map is 1:150,000,000. I imported the Thunderforest Pioneer base layer. At the scale I used, country labels do not show when working on the map in Print Composer.

When I export the map as a PDF, png or svg, the country labels now appear. I do not want them to. How can I export the map as is being displayed in Print Composer. The basemap does not have an option to modify the labels.

Here's a screen shot of what I see in Print Composer:

Note that I locked the layers and locked the layer style.

Here's a screen shot of what I see in the PDF export:

Here's a PDF of what I see when I right click on the base layer and select properties:


texts: How can I get the plugin manager "get more" icon?

texts:
I am fairly new to QGIS and i was looking to add more plugins to the core provided in version 2.0
I read the user manual provided by QGIS and found out that there should be an icon "get more" available in plugin manager. Actually in the QGIS version i am using (2.0.1 Dufour) I can only see two icons: installed and settings.
How can i get the get more one? Any suggestion?
In addition, if i go in settings the status of the plugin repositories is unavailable. Is this "normal"?

I have checked no Get More in Plugin Manager in qgis 2 + ubuntu 13.04 but I get a diffeterent window for the settings.


texts: How to save point ID when extracting values from bands?

texts:
I have vector points that have some external IDs. I would like the point ID column to appear in the result table after extracting spectral values from the bands. 'Latitude', 'longitude', 'time' and '.geo' columns can be skipped, but I see that they are mandatory from the getRegion() function.

var points = [ee.Geometry.Point(-77.64, 37.51),
              ee.Geometry.Point(-77.66, 37.53),
              ee.Geometry.Point(-77.65, 37.52),
              ee.Geometry.Point(-77.68, 37.52),
              ee.Geometry.Point(-77.62, 37.54)];
var table = ee.FeatureCollection(points);

// external points ID
var dict = {oid:[101, 102, 103, 104, 105]};
var table = table.set(dict);

var data = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
              .filterDate('2018-01-01', '2018-04-01')
              .filterBounds(table)
              .getRegion(table, 30);

var keys = ee.List(data.get(0));
var featureCollection = ee.FeatureCollection(data.slice(1).map(function(singleData){
  singleData = ee.List(singleData);
  var dict = ee.Dictionary.fromLists(keys, singleData);
  var point = ee.Geometry.Point([dict.get('longitude'), dict.get('latitude')]);
  return ee.Feature(point, dict);
}));

Export.table.toDrive({
  collection: featureCollection,
  description: 'CSV'
});

Link to code


texts: Earth Engine sample regions functions after a loop returns different formatting than individual output within the loop

texts:
I have a feature collection of points which has a date column within it.  I want to loop through the unique dates in the feature collection and then match corresponding dates within an image collection.  After matching the features and images I want to extract values from the images using the sampleRegions function.

I have successfully done this manually with one date at the top of my script, and then tried to implement the same thing in a loop which seems to return seemingly empty output.

Here is a link to my code:
https://code.earthengine.google.com/ea5f56d6c7b98619fbeffa373c9d95ed

And here is the code:

var bui = ee.ImageCollection("users/spotter/GFWED_v25/MERRA2_BUI"),
    pts = ee.FeatureCollection("users/spotter/ak_ca_combustion/raw/sol_pts");

//print a single pts feature to see what the first Date column item is
print('First Feature', pts.first()); //date is '2001-04-24'

//------------------------this works for extracting a single date without a loop
//the same date in the first feature
var date = ee.Date('2001-04-24');

//filter the pts by this date
var sub_pts = pts.filter(ee.Filter.eq("Date", '2001-04-24'));

//filter the image collection by this date
var bui_sub = bui.filterDate(ee.Date(date), ee.Date(date).advance(1,'day')).mean();

//extract the pts that match the image collection by date
var extract_test = bui_sub.sampleRegions({collection: sub_pts,
                                scale: 460,
                                projection: 'SR-ORG:6974',
                                tileScale: 16,
                                geometries: false
  });

print('Extracted Data', extract_test);

//---------Now I want to do the same thing over many dates at onc

//select specific dates of interest
var all_dates = ['2001-04-27', '2001-04-28', '2001-04-29'];

//loop many dates at once

var feat_iteration = ee.List(all_dates.map(function(date){

  //get sub image collection matching date
  var bui_sub = bui.filterDate(ee.Date(date), ee.Date(date).advance(1,'day')).mean();

  //get sub pts matching date
  var sub_pts = pts.filter(ee.Filter.eq("Date", date));

  //get the values for each point feature in sub_pts to be extracted from bui_sub
  var extract = bui_sub.sampleRegions({collection: sub_pts,
                                scale: 460,
                                projection: 'SR-ORG:6974',
                                tileScale: 16,
                                geometries: false
  });

  //this looks fine
  print('Each Extract', extract);
  return(extract);
}));

//the output data here is formatted much differen't than individual print statements
print(ee.FeatureCollection(feat_iteration));


texts: How can I linearly order KML placemarks?

texts:
Is there anyway to take an unordered set of placemarks in a KML file and add a linear order to them?

By linear order I mean add a numeric number so that each placemark has an order in the group of placemarks.  This will essentially turn the placemarks into a ordered path.

I'm willing to use any program to import the KML file to do these operations however it has to have a GUI so that semi-technical people can do it. (i.e. people without knowledge of xml)


texts: Setting up ArcPad applet (Python) that calls arcgisscripting?

texts:
I've been asked to support a legacy script that was "made" (read: found) by a former coworker.

My python skills are pretty basic, so I went right to installing the applet like I would for a VBS applet, got an error saying no python, so I put on python, similar error, I read about python win32, and got that on all nicely matched up with bit (32) and version (2.7.7).  Ran the applet, and it starts python and does some string manipulations to create an output name so I know Python is working, then it crashes on "import arcgisscripting".

Now I have enough Python experience to know that's the usual call you make in the Desktop environment... A quick check and sure enough the laptop where this applet currently works has a copy of ArcGIS on it...

I have read in a few places that a ArcGIS runtime SDK might solve the problem, but which one...?  Does python even care?  How do you buy a runtime sdk license and how much?

Many many thanks for any help at all.  At this point I'm open to any suggestions.


texts: Is there an Offset function in QGIS for Polygons?

texts:
IS there an offset function in qgis similar to that in autocad for polygons. I have an irregular shaped pollygon which i want to offset at 100m out from the edges, effectively adding 100m out each time.

I can do it buffering but does not show the exact shape.


texts: Filtering a Feature Class in using ArcGISDynamicMapServiceLayer

texts:
I'm using the following code in Flex to display all the features in a Feature Class defined in the "Data Layer".

var dataLayer:ArcGISDynamicMapServiceLayer= map.getLayer("Data Layers") as ArcGISDynamicMapServiceLayer;
var dataLayerId:int;

dataLayerId = 5;

dataLayer.visibleLayers.removeAll();
dataLayer.visibleLayers.addItem(dataLayerId);
dataLayer.refresh();

But I'd like to apply a filter like "WHERE MGD = 1" (MGD is an attribute defined in the feature class) to display only a subset of that feature class.
Any help or autosuggestion would be much appreciated.


texts: BMP Siting Tool gives error at validateDatasets_ST: Error HRESULT E_FAIL returned from call to COM component?

texts:
I'm trying to use a toolbar that was developed by the EPA that is BMP Siting Tool. I also have done what was ordered in BMP Siting Tool User Guide, but still failed until now.
When I run BMP Siting Tool and wait for the result, then a dialog box appears

  error at validateDatasets_ST: Error HRESULT E_FAIL has been returned from a call to a COM component

Any idea what's causing this? I'm using ArcGIS 10.1. If maybe anyone have a sample file for using on BMP Siting Tool?


texts: histograms chart widget in CMV application

texts: I am trying to build my first widget te creat charts or histograms, for use in the Configurable Map Viewer (CMV) im asking if you can help me with tutorial or something else can help me to build this  widget in the CMV application


texts: I have shapefile with roads. I want to digitize some of the objects in shapefile to OSM, because I allways use OSM shapefiles and I have something to give to. What software should I use? I have JOSM installed. Is it enough?


texts: QgsVariableEditorWidget does not save changes

texts:
I would like to edit the variables of a composer item in QGIS 2.18.x.
The widget itself is working so far. In my example, I am selecting a QgsComposerLabel called "Plankopf" to show its variables.

What is not working so far is:

saving the edits: I thought that changes done throught the QgsVariableEditorWidget would be saved but that's not the case.
What do I have to change to get added/updated variables saved?
showing changed Variables in the Composer-Window: Even if I change variables with QgsExpressionContextUtils.setComposerItemVariable(), the changes are not visible in the Composer window. How can I refresh the variables widget in the Composer window?

# -*- coding:utf-8 -*-
from qgis.core import QgsExpressionContextUtils, QgsExpressionContext
from qgis.gui import QgsVariableEditorWidget
from qgis.utils import iface
from PyQt4.QtGui import QInputDialog, QWidget, QFormLayout, QPushButton, QLineEdit
from PyQt4 import QtCore
from PyQt4.QtCore import QPoint

class inputdialogdemo(QWidget):
    def __init__(self, parent = None):
       super(inputdialogdemo, self).__init__(parent)

    def getItem(self):
       items = []
       composerId = 0
       composers = iface.activeComposers()
       for item in composers:
           items.append( item.composerWindow().windowTitle() )
       item, ok = QInputDialog.getItem(self, u'Auswahl',
         u'Druckzusammenstellung wählen', items, 0, False)
       if ok and item:
           return item

 ex = inputdialogdemo().getItem()
 print ex

 composerId = 0
 composers = iface.activeComposers()
 for item in composers:
     if item.composerWindow().windowTitle() == ex:
         break
     composerId += 1

target_citem = NULL
def itemsuche(name):
    for citems in composers[composerId].composition().items():
        if type(citems).__name__=='QgsComposerLabel':
            if citems.displayName() == name:
                return citems
                break

target = itemsuche('Plankopf')

canvas = iface.mapCanvas()
variable_editor_widget = QgsVariableEditorWidget()
expression_context = QgsExpressionContext()

expression_context.appendScope(QgsExpressionContextUtils.composerItemScope(target))

variable_editor_widget.setContext(expression_context)

variable_editor_widget.reloadContext()
variable_editor_widget.setEditableScopeIndex(0)
print(variable_editor_widget.context())
print(variable_editor_widget.editableScope())
print(variable_editor_widget.settingGroup())
print(variable_editor_widget.variablesInActiveScope())

mapcanvas_upperleft = canvas.mapToGlobal(QPoint(0, 0))

offset_x = 0
offset_y = 30
pos_x = mapcanvas_upperleft.x() + offset_x

pos_y = mapcanvas_upperleft.y() + offset_y
variable_editor_widget.setGeometry(pos_x, pos_y, canvas.width()/3, canvas.height()*0.75)

variable_editor_widget.show()
print variable_editor_widget.pos().x()
print variable_editor_widget.pos().y()
target.setHtmlState(0)
target.setHtmlState(2)


texts: Changing Opacity of Vector layer in Google Earth Engine

texts:
While importing a Shp file in Earth Engine, I want to control or change opacity/transparency of the layer. I have gone through almost all read me files, I find nothing.

This doesn't work:

Map.addLayer (ind_shp, {opacity: 0.5}, 'map');


texts:
I have a table of Census Blocks linked to ZIP4 values. I'm trying to create another table getting the unique ZIP4 values associated with each Census Block.

I've tried using the Frequency, and Summary Statistics tool. These give me what I need, but the ZIP4 values are all added up, I need them separately.

Would using Pivot Tables get me what I need?

Using ArcMap 10.2.2, all extensions.


texts: Move features to different geographic area

texts:
I want to move a dataset of a few thousand vectors to a different geographic area. However I can't seem to find any tools to do this, or I'm not sure on the terminology for what to look up.

I suppose a common example of what I would like are maps that show the true size of Antarctica when overlaying other continents. However, size comparison isn't my actual goal if it makes a difference; any size, direction, or shape distortion isn't a concern. My intent is to obfuscate the real location while still being able to show geographic features.

My preferred tool is QGIS but I can use ArcMap or PostGIS.


texts: Analyzing signatures for Landsat clasification with Sextante in QGIS

texts:
I am defining my training data for a Landsat classification. I am digitizing the polygons with the field data in QGIS and I want to analyze the signatures of the Landsat images to avoid errors.
In ERDAS I use the Signature editor, but I think I can do it with the Sextante toolbox in QGIS.
Does anyone know of an example of this on the web?


texts: Doing NDVI when nir data is given as 3 bands

texts:
I have been provided with two data sets of the same area, one with normal RGB bands and one with NIR data. My problem is that since the NIR data is given as 3 bands, i am not sure how I should interpret this.

I am guessing that its normal in some data acquisition areas to collect NIR data with a DSL camera with some filters removed and thats why the data is captured over 3 bands. (but i really dont know much about this).

How do i transform the 3 bands of nir data into the one I need for my NDVI calculation?


texts: How do we use the agr argument in st_sf()?

texts:
I am using the sf package to create new spatial datasets, via the st_sf() function. I do not understand the use of the attribute-geometry relationship argument of this function. The help file specifies that it defines how each non-geometry field relates to the geometry field, and takes one of three values for each row: "constant", "aggregate", or "identity". There is also a value for unknown.
Here is what I can understand from the help file:

"constant" is used for attributes that are constant throughout the geometry (e.g. land use)

For a polygon or a line, this indicates that whatever the attribute is describing does not change inside the feature, such as (per the example) a homogeneous polygon of land use.

"aggregate" where the attribute is an aggregate value over the geometry (e.g. population density or population count)

From what I can tell from the example, this is used when the attribute is calculated from the geometry itself. Area and perimeter in the sf example are classified as 'aggregate' attributes.

"identity" when the attributes uniquely identifies the geometry of particular "thing", such as a building ID or a city name

And with respect to "identity", I cannot find a salient example of this being different from "constant".
To test this, I loaded demo(nc) and tried various combinations of agr values. Changing the values randomly had no impact on the class of the fields, their values, or how they behaved when summarized.
Can someone please explain why this argument exists, and how (or if) I should be using it when creating features?


texts: Problem with input feature as clipping geometry in ArcGIS Pro

texts:
I am clipping several rasters to the shape of a polygon feature class using the feature class for clipping geometry. However the output looks like this:

There should be NoData values around the polygon instead of zeros. Any idea about this problem?

I am using ArcGIS Pro 2.4.2


texts: Volume calculation of different termite mouts

texts:
So i know how to calculate the volume of 1 indiviual mountain, but my problem is that i have a data set (LiDAR) with a lot of those. I was thinking if there was a more effecient way to calculate the volumes.
I made a pointcloud in Agisoft Metashape, then moved it to cloudcompare as a .las file, then made a DEM and DTM and moved it to QGIS and calculated a nDSM with contours.
2 things i think of if possible:

Is it possible to extract the mounts based on the height (i made the DTM so every mount has a comparable starting hight.
Can you calculate the mounts separate (if i try it now it just gives me 1 total volume)

If anyone thinks of a better way let me know.


texts: Setting the CRS of Raster in standalone PyGQIS

texts:
I'm generating raster files without a CRS and would like to assign one. I am attempting to use GDAL's assign projection, which works perfectly inside QGIS, but seems to have no effect when run in my standalone python environment. Any thoughts?
    alg_params_DEM = {
    'ADVANCED_MODIFIERS': '',
    'CELLSIZE': 1,
    'CLASS': '2',
    'INPUT': las_loc,
    'MEDIAN': '',
    'MINIMUM': '',
    'SLOPE': '',
    'SMOOTH': '',
    'SPIKE': '',
    'VERSION64': True,
    'XYUNITS': 0,
    'ZUNITS': 0,
    'OUTPUT_DTM': home_folder+'\\DEM.dtm'
    }

    dem = processing.run('fusion:gridsurfacecreate', alg_params_DEM)

    processing.run("fusion:dtm2tif",
    {'INPUT':home_folder+'/DEM.dtm','MASK':False,'OUTPUT':home_folder+'/DEM.tif'})

    processing.run("gdal:assignprojection",
    {'INPUT':home_folder+'/DEM.tif','CRS':'QgsCoordinateReferenceSystem("EPSG:28355")'})


texts: Determining the transformation between a satellite image and long/lat

texts:
I'm trying to overlay some data onto a satellite image of the UK.

The image is a hi-resolution version of the image below (originally from dundee university):
http://www.bbc.co.uk/news/uk-scotland-18239311

However, I'm unable to determine the mapping between longitude/latitude and points on the image, i.e. given a long/lat pair what is the relevant point on the image. Along the way I created a table of landmarks round the edge of the UK and recorded their position on the image (with top-left being (0,0)) together with their longitude and latitude values:
https://gist.github.com/pci/4945162 . Sadly I don't know the location of the satellite or the direction the camera was pointing when the shot was taken.

I've tried affine transforms, perspective transforms and planar projections from a sphere, all with very little success.

Any help is greatly appreciated.


texts: How to let rasterio use a colormap file when showing a TIFF

texts:
I'm new to rasterio. I have tried this example. My raster file is a GeoTIFF file and I have a colormap file available in ESRI *.clr format, with 58 lines. The values in the GeoTIFF file range from 10471 thru 10534. For those who don't know the format, the first 2 lines looks like this:
10471 16 24 79
10472 104 223 247
I have read a bit about colormaps, e.g. here. I read this: "You cannot set the color interpretation on existing data". All the same, I'm hoping that I can have the colormap applied to the plotted raster in a direct way, but I have no idea how to go about it. In order to start somewhere, I wrote this code:

import csv
import os.path
import sys
cf = None
try:
    curdir = os.path.dirname(sys.argv[0])
    datapath = os.path.normpath(os.path.join(curdir, "../geodata/EUR"))
    cf = open(os.path.join(datapath, "countries.clr"), "r")
    rows = {}
    for row in (csv.reader(cf, delimiter=' ')):
        color = tuple([int(row[1]), int(row[2]), int(row[3]), 255])
        rows[int(row[0])] = color
    # Do we almost have a colormap now?
finally:
    if cf != None: cf.close()

Later on, I'd also like to indicate that rasterio should plot classes - defined somehow - in different colours. I'm afraid that it's only possible to get this working in an indirect way - with some time-consuming tricks - but if you can, please surprise me!


texts: How to compute flow direction and flow accumulation with Qgis?

texts: Any one who had made flow direction and flow accumulation on raster data for filled evaluation model i Qgis? Which tool and plugins are necessary, i haven't been able to find any guidelines on Google


texts: Deleting features (rows) from shapefile (attribute table) using ModelBuilder?

texts:
I've run into a problem and can't figure out how to solve this one.

I used the ModelBuilder and have a shapefile with loads of information. This information needs to be reduced. Therefore I want to select several rows (not all of them) out of several columns by attribute to delete them.

So for example, column X contains data with attributes a, b and c. I only want to delete c. Since I'm using the ModelBuilder, I need a tool from the toolbox which does that (preferably) in one step.

Any ideas?


texts: coordinates point in any reference system

texts:
I have a coordinates point like this "01010000003413056727B71D41C7F1A68B1ADC4141", but I don't know it is what.

Can you aswer What is this?


texts: Units of curvature output in ArcGIS 10.x and Pro

texts:
This question has been asked several times, but I am not sure that there is a satisfactory answer.

I have used ArcGIS Pro to calculate curvature from a DEM. Now I am not sure what the units are in the curvature raster. According to ArcGIS, the curvature values are calculated using the second derivatives of the input surface in a 3x3 window. Here is an example that comes from my dem. The average elevation is roughly 325 meters, with a cell size of 1 meter. The curvature value for the center cell is -28.8. I obtain the same output using the following code in python. My question is, what does -28.8 (or -0.288, if you don't multiply by 100) mean?

df = pd.DataFrame(np.array([[325.88,325.66,325.73],[325.799,325.670,325.700],[325.920,325.809,325.850]]))

Z1 = df.iloc[0,0]
Z2 = df.iloc[0,1]
Z3 = df.iloc[0,2]
Z4 = df.iloc[1,0]
Z5 = df.iloc[1,1]
Z6 = df.iloc[1,2]
Z7 = df.iloc[2,0]
Z8 = df.iloc[2,1]
Z9 = df.iloc[2,2]

A = ((Z1 + Z3 + Z7 + Z9)/2 - (Z2 + Z4 + Z6 + Z8)/2 + Z5) / 1**4
B = ((Z1 + Z3 - Z7 - Z9)/4 - (Z2 - Z8)/2)/ 1**3
C = ((-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6) /2) / 1**3
D = ((Z4+Z6)/2 - Z5) / 1**2
E = ((Z2+Z8)/2 - Z5) / 1**2
F = (-Z1 + Z3 + Z7 - Z9) / (4*1**2)
G = (-Z4 + Z6) / 2*1

curvature = -2 *(D + E)*100

If you change the values to df = pd.DataFrame(np.array([[1,1,1],[2,2,2],[3,3,3]])), curvature is 0, as expected. If you change the values to

df = pd.DataFrame(np.array([[1,1,1],[2,2,2],[4,4,4]]))

, the value for curvature becomes 100. Again, 100 (or 1) what?


texts: Burning stream network into DEM layer using ArcGIS Desktop?

texts:
I am working on a DEM that in an area where there is little relief.

I would like to 'burn' the river network into the DEM so I can calculate flow accumulation and flow length accurately.

I am using ArcGIS Desktop 10.


texts: Novel ways of filling in missing values in a raster

texts:
I am looking for robust ways to fill in missing values in some rasters. They all have a single layer. Missing values consist of single pixels to medium sized patches. Rasters are around 1000 x 1000 pixels in size and the largest patches are like 20x20 pixels.

I'm tempted to use aregImpute in the Hmisc R package. Has anyone used it for this purpose?

this approach looks seems very cool but I think it is only meant to produce aesthetically pleasing corrections.

Detailed explanation of this:

All the rasters (I have 36 in total) share the same extent, they overlap and are aligned. Each raster is a different variable, I gathered variable from various sources (remote sensing, topographic and climatological). The original rasters come in various resolutions. The smallest being 30m. From there up they get as high as 1km. I resampled everything using cubic convolution (all the variables are continuos) to 1km. I have another 1km raster where I have data of a variable of interest for some sampled points. So I trained a model using those points and the other rasters as covariates to be able to generate a full raster of that variable. Unfortunately most covariate rasters have some missing values in them, actually not much but I would want to eliminate the problem entirely.

Thank you.

p.s. I would preferably like to use R for this.


texts:
I have several (100+) .kmz files that I would like to import into my PostgreSQL database. I have this ogr2ogr command that works fine getting it into the database, but creates a new table every time:

C:\Radio Stations>ogr2ogr -f "PostgreSQL" PG:"host=localhost user=postgres dbname=postgres password=password port=5433" "file.kmz"

What do I put in that cmd to put all of the .kmz files into the same table?

I'm using PostGIS 2.2.1, PostgreSQL 9.5.1, Python 3.6.1


texts: Rotating feature perpendicular to feature in another layer using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop basic licence I have created polygons (rectangles) in one layer and i want to move them so that they are exactly perpendicular to a boundary in another layer.

I have been using the rotate tool and 'a' to enter specific angles but it is quite time consuming. I know I could have created the polys using the sketch tool and aligning it perpendicular to the boundary layer.

but i wanted them exactly the same size, distance, orientation etc to each other, so i simply copied and pasted.
I am trying to avoid restarting. There is quite a lot on here r.e. rotation but not exactly what i need.


texts:
I need to get data for my database from some GRIB files (binary), I downloaded wgrib2, which can export GRIB-2 files to .csv and csv is very suitable for me, but some GRIB files are GRIB-1, so I used wgrib2 on them and program wgrib2 told me: "grib1 message ignored (use wgrib)", so I used wgrib, but wgrib can't export files to different format.
So I tried to install recommended cnvgrib to convert GRIB-1 to GRIB-2, but still getting some errors during installation due different library versions, 32/64 bit problem, different directories etc., it's written in Fortran, that's also the problem I can't compile it easily...

So is there another way to convert GRIB-1 files to GRIB-2 (or to some non-binary format)? I can google many ways to do it, but I need help from some1 who have already converted GRIB-1 to another format, because there can be some easy way how to do that, but I have so much pain fixing installation problems, wrong configurations, Fortran version incompatibilities etc during installation of cvngrib, which I haven't even fixed and gave up...

EDIT:
OK, my co-employee fixed degrib (it linked own libraries with newer system libraries), so I can run degrib now.
So I run "degrib "GRIB file" -C -msg all -nMet -Csv" on GRIB-1 file and it failed the conversion to .csv with:

ERROR: In call to Grib2Convert.
ERROR: In call to ReadGrib2Record.
Problems with ReadGrib1Record called by ReadGrib2Record
Inside ReadGrib1Record
Grid projection number is 4
Don't know how to handle this grid projection.

I try to fix it soon, but meanwhile you can help me again!


texts: Clipping object of class LAS with polygon shapefile using lidR?

texts:
I want to make a simple clip with my .las 3D-Pointcloud.

Therefor I use the lasclip() function.

Import .las:

lidar <- readLAS("~/data.las")
lidar

gives me

class        : LAS (LASF v1.2)
memory       : 774.7 Mb
extent       : 3403298 , 3403640 , 5286553 , 5286831 (xmin, xmax, ymin, ymax)
area         : 62291.4 m² (convex hull)
points       : 11946448 points,  4704414 pulses
density      : 191.78 points/m²,  75.52 pulses/m²
field names  : X Y Z gpstime Intensity ReturnNumber NumberOfReturns Classification ScanAngle R G B pulseID
coord. ref.  : +init=epsg:31467 +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7

Import the clip:

site <- readOGR(dsn = "~/directory", layer = "clipshape")
site

gives me:

class       : SpatialPolygonsDataFrame
features    : 1
extent      : 3403372, 3403590, 5286594, 5286817  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
variables   : 1
names       : id
value       :  1

Ok now I want to clip it:

lasclip(lidar, site)

and I get the following Error:

Error: Geometry not supported

I am working in an environment where I cannot use pdal with python for several reasons.

How can I get cliplas() working?


texts: Automatically display related records for selected feature in tablewindow with ArcObjects

texts:
Is there a way through ArcObjects to automatically display the related table(s) for the selected feature in ArcMap in the attribute window (tablewindow) without the user having to click on it and select the related table. I can add a button that displays the tablewindow but I don't see how I can automatically display the related tables tab that appear in the window when the "related tables" button is clicked.

There is an ITableViewCallback interface that has a ShowRelationshipTable method to display related records. I can create a class that implements that interface but that doesn't seem to be what I want. I really just want to call the "ShowRelationshipTable" method of the tablewindow which doesn't seem to be exposed. I've also tried executing the command for "Related Tables" in the tablewindow which I think is:
esriGeoDatabaseUI.TableContextMenuRelationships but that doesn't do anything either. Is there a way with ArcObjects to automatically display the tabs for the related data without the user clicking on it?


texts: Does PostGIS Catch Self-Intersections Involving the Start / End points to a Line or Area?

texts:
I'm trying to see if I can catch among other things a self-intersection which involves a start / end point of a geometry feature self-intersecting some other part of the same feature.  Obviously the start and end points of a proper feature will be the same.

I ask because I've heard some tools ignore start and end points altogether when trying to determine a self-intersection.  This can be a problem if either start or end points legitimately intersect some other point / line segment.

Would catching this specific kind of self-intersection be easily caught via ST_IsValid?


texts: Advanced branching in Survey123 for ArcGIS

texts:
I'm designing a webform in Survey123 for ArcGIS which requires those responding to the survey to identify bus stops along a route, one by one, to essentially submit the full sequence of stops so that the bus route can then be mapped.

To simplify, the flow of the survey is as follows:

The respondent selects a bus service from a drop-down of pre-defined options (the list of selectable services will be derived from a previous survey delivered through Survey Monkey).
The respondent then identifies which stop along the route sequence they are about to input. The stop number is selected from a drop-down list of pre-defined responses.
Finally, the respondent then provides the stop location by entering a data point on the Esri-powered geopoint map widget.
After submitting the data point, the respondent has the option to return to the form and select the next stop along the route and then once again provides the location on the geopoint widget. The idea is that the user will keep returning to the form to input data points until they have entered all points along a bus service.

Ideally, there would be a singleline-text question between 1 and 2 asking the respondent how many stops there are along their chosen service from Step 1. However, there does not seem to be a way to link such a singleline-text question so that the numerical value inputted into it would inform the possible
options in the drop-down question at Step 2. As such, i have had to set up Step 2 so that the pre-defined options range from 1-30 as i do not know the highest number of stops amongst the bus services, and a maximum of 30 should cover all possible ranges of numbers of stops as the services are expected to be quite rural (although i cannot be sure at this stage where the services go through, nor how frequent stops are along them).

I have been informed that setting up a drop-down question to have options which are defined by setting a conditional relationship with a previous singleline-text question is possible within Survey Monkey, and that this
relationship is referred to as "advanced branching," which from my understanding dictates how often a respondent gets asked a question depending on a value they provided earlier in the survey. I've also been told that it is not possible to build advanced branching into a Survey123 form as the feature class in a geodatabase which has been generated for Survey123 must have a set number of fields defined before data entries are added to it. Is this true? Can i design the survey so that the respondent can only select stop numbers (and assign locations for said stop numbers) based
on the total amount of stops they will have identified in the previous answer within the survey?

While i have been using the web browser version of Survey123 for ArcGIS to design the webform, if this functionality can only be added through the advanced desktop version (Survey123 Connect) then i have been trying it out and have some proficiency with that version.


texts: I've mapped about 10,000 pick up and the corresponding drop off points. I wanted to connect each pick up point to it's respective drop off point with a line. What is the best way to go about doing this in QGIS? If this is not practical or if there is a better visualization method that someone knows of I am welcome to recommendations as well. I have tried using the plugin MMQGIS and its function Create - Hub Lines. I have also tried using Points2One function but none of them seem to do what I want (which I know from google searching has been done before, although on a smaller scale).


texts:
I would like to use QGIS with QtQuick. I can use QGIS with QtWidget with the libs that I downloaded from OSGeo4W but there is no qtquick libraries in there. In QGIS's website there is a demo with QtQuick but since I don't have the qtquick libraries I can't use that code. Where can I found qgis quick libraries?

In demo's profile they write this but I don't have qgis_guick lib

  QGIS_PREFIX_PATH = $${QGIS_INSTALL_PATH}
  QGIS_LIB_DIR = $${QGIS_INSTALL_PATH}/lib
  QGIS_INCLUDE_DIR = $${QGIS_INSTALL_PATH}/include
  QGIS_QML_DIR = $${QGIS_INSTALL_PATH}/qml
  QT_LIBS_DIR = $$dirname(QMAKE_QMAKE)/../lib

  INCLUDEPATH += \
    $${QGIS_INCLUDE_DIR}/qgis \
    $${QGIS_INCLUDE_DIR}

  LIBS += -L$${QGIS_LIB_DIR}
  LIBS += -lqgis_core -lqgis_quick


texts: Exporting tiff for use in photoshop problems

texts:
I am working with a section of the SRTM DEM in .tiff format (8bit unsigned integer) in ArcGIS 10.3. I want to export the data in order to do some modifications in photoshop. However when I export data (tried different formats and compressions) and try to open it in photoshop I keep getting completely black or a partially black images in photoshop.

Does someone know where I am doing something wrong?


texts: Installing GDAL on Windows 10 using binaries for command line functionality

texts:
I figured that I would ask this question hoping that someone finds a document or provides a description that actually works. I've followed countless sites, GitHub issues, and other forum posts but have not been able to preserve command line functionality with Python scripts.
Using: Windows 10 and Python 3.7, 64 bit
I do not want to use OSGeo4W or their shell, or conda. I need to use the Windows binaries. I also need to use GDAL version 3.1 due to options not found in GDAL version 3.0
Here is what I have attempted so far. This installation is from scratch, without any dependencies installed.

Navigate to the binaries (this version) and download the following: gdal-300-1911-x64-core.msi for the core (do typical install) and GDAL-3.1.0.win-amd64-py3.7.msi for the python bindings.
Add C:\Program Files\GDAL to your PATH under system variables. Move it to the very top.
Add these two system variables:
GDAL_DRIVER_PATH=C:\Program Files\GDAL\gdalplugins and GDAL_DATA=C:\Program Files\GDAL\gdal-data
Press OK and navigate to your command prompt.
Type in gdalinfo --version. This should output GDAL 3.1.0dev, released 2019/99/99
Type in gdal2tiles.py --xyz

At step 6, the application will break. Obviously, we need an input file, but that doesn't matter right now because we receive a different error. Here is what the command prompt outputs:
C:\Users\Bold>gdalinfo --version
GDAL 3.1.0dev, released 2019/99/99

C:\Users\Bold>gdal2tiles.py --xyz
Traceback (most recent call last):
  File "C:\Users\Bold\AppData\Local\Programs\Python\Python37\lib\site-packages\osgeo\__init__.py", line 18, in swig_import_helper
    return importlib.import_module(mname)
  File "C:\Users\Bold\AppData\Local\Programs\Python\Python37\lib\importlib\__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "<frozen importlib._bootstrap>", line 1006, in _gcd_import
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 670, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 583, in module_from_spec
  File "<frozen importlib._bootstrap_external>", line 1043, in create_module
  File "<frozen importlib._bootstrap>", line 219, in _call_with_frames_removed
ImportError: DLL load failed: The specified module could not be found.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Program Files\GDAL\gdal2tiles.py", line 53, in <module>
    from osgeo import gdal
  File "C:\Users\Bold\AppData\Local\Programs\Python\Python37\lib\site-packages\osgeo\__init__.py", line 32, in <module>
    _gdal = swig_import_helper()
  File "C:\Users\Bold\AppData\Local\Programs\Python\Python37\lib\site-packages\osgeo\__init__.py", line 31, in swig_import_helper
    return importlib.import_module('_gdal')
  File "C:\Users\Bold\AppData\Local\Programs\Python\Python37\lib\importlib\__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
ImportError: DLL load failed: The specified module could not be found.

As you can see, we receive a weird error. The .exe files work, but not the .py ones. It's most likely due to the Python bindings, but I believe that they are set up correctly.


texts: How to merge layers from two vector tiles with Mapnik?

texts:
I'm using Mapnik + Mapbox nodejs stack.  Given two vector tiles (PBFs) for the same location/zoom, with some layers having the same name in both. How do I merge them into one PBF tile - replacing identically named layers in the first tile with the ones from the second tile?

Update:
I found that I can add layers via GeoJson, but it seems awfully inefficient: addGeoJSON does not replace existing layer, but adds a new one (even if the same name already exists), so I have to export each layer to JSON and re-add them to a new tile:

newVectorTile.addGeoJSON(vectorTile.toGeoJSON('water'), 'water');


texts: Open the popups of all markers that are visible on the map with a zoom > 10

texts:
I need to open all the markers that are visible on the map at Zoom > 10. I am using leaflet.markercluster plugin.

Init map:

initMap() {
  this.map = L.map('map', {
    center: [55.55, 37.61],
    zoom: 9,
    layers: this.layer
  })
  this.tileLayer = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution:
      '&copy; OpenStreetMap, &copy;'
  })
  this.tileLayer.addTo(this.map)

this.map.on('zoom', function(ev) {
    ???
  })

Add marker layer:

this.markerLayer = new L.LayerGroup()   // layer contain searched elements
  // this.map.addLayer(this.markerLayer)

  for (const i in data) {
...
    const marker = new L.Marker(new L.latLng(loc), { title: title, icon: icon })// se property searched
    marker.bindPopup(title)
    this.markerLayer.addLayer(marker)
  }

Use leaflet marker cluster:

this.markersLayer = L.markerClusterGroup({
    iconCreateFunction: function(cluster) { ... },
    singleMarkerMode: false
  })
  this.markersLayer.addLayer(this.markerLayer)
  this.map.addLayer(this.markersLayer)


texts:
So, I have a polygon that contains the population number. Another polygon is a zone polygon.
I need to approximate population from the first polygon to zones polygon, but based on intersection percent.

Here is an example. Purple polygon contains population data. Now, how to split for example population of 64 to other polygons based on the percent of their intersection.
I tried to use ratio policy on a polygon with population data, then clip it with other polygon, but with no luck.


texts: How to convert SPOT NDVI hdf file to tiff file using gdal_translate command

texts:
I have SPOT NDVI images. I want to convert NDVI hdf to tiff. I have used the following code to get the information of hdf file:

gdalinfo 20060501_NDVI.HDF

I get the following output

Driver: HDF4Image/HDF4 Dataset
Files: 20060501_NDVI.HDF
Size is 8849, 5601
Coordinate System is `'
Metadata:
 valid_range=0, 255
Corner Coordinates:
  Upper Left  (    0.0,    0.0)
  Lower Left  (    0.0, 5601.0)
  Upper Right ( 8849.0,    0.0)
  Lower Right ( 8849.0, 5601.0)
  Center      ( 4424.5, 2800.5)
  Band 1 Block=8849x113 Type=Byte, ColorInterp=Gray

As I do not get any information about the SUBDATA SETS in the gdalinfo command I use the following code to convert HDF to TIFF

gdal_translate -of GTiff 'HDF4Image:"20060501_NDVI.HDF":*NDVI'

But I get the output as

FAILURE: No target dataset specified.

As no DATASET is specified, I unable to specify any target dataset. Please let me know how to resolve this. Also I would like to resample the image to 250m resolution by using nearest neighbour (default) method.


texts:
Is it possible to determine whether a point feature of a particular height would be backed by sky or by terrain?

I've done the viewsheed analysis of a wind turbine. Now I would like to make a raster map representing the background (sky/terrain) from each cell of the viewshed.


texts: I am new to QGIS 1.8 and have been trying to make a heatmap.  From what I can tell, the reason I am not getting the results I want is that the max value of the raster is erroneously measured at 2.  Max should be 15-20 from what I can tell. I have tried to run the operation multiple times with different parameters and the max value is 2 every time. I have used the Identify tool to verify that the problem is not the cells themselves. Any cell with value above 2 is symbolized the same. Any ideas on how to fix this?


texts: Converting Digital Number (DN) to reflectance in order to analyze NDVI

texts:
I'm working with a mapir nir + red camera meaning that its blue band captures nir light and its red captures red normally. I understand that NDVI is calculated using reflectance values.
How do I convert 8 bit RGB into reflectance values so that I'm able to calculate each NDVI value? Also, does converting an image to grayscale mean its reflectance?


texts: Debugging FeatureClassToFeatureClass giving ERROR 000732?

texts:
I get this Error when trying to run the FeatureClassToFeatureClass_conversion:

  Traceback (most recent call last):   File "", line 1, in

      arcpy.FeatureClassToFeatureClass_conversion(inFeatures, outLocation, outFeatureclass)   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\conversion.py", line 1675, in
  FeatureClassToFeatureClass
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Features: Dataset Parcel_Owners does not exist or
  is not supported Failed to execute (FeatureClassToFeatureClass).

This is my code:

import arcpy
arcpy.env.workspace = "D:\Users\ShawnO\Desktop\PythonExercizes\Parcel_Ownership_update_model\Parcels.gdb"
inFeatures = "D:\Users\ShawnO\Desktop\PythonExercizes\Parcel_Ownership_update_model\Parcels.gdb\Parcel_Owners"
outLocation = "D:\Users\ShawnO\Desktop\PythonExercizes\Parcel_Ownership_update_model\Old_Parcel_Owner\Old_Parcel_Owners.gdb"
outFeatureclass = "old_owners"
arcpy.FeatureClassToFeatureClass_conversion(inFeatures, outLocation, outFeatureclass)

What is wrong with my variable inFeatures? I know that a feature class table consists of more than one file but examples on the internet had the name of the Featureclass as an acceptable input.


texts: QGIS Raster merging UnicodeDecodeError: invalid start byte

texts:
Edit: I moved them to another folder so that the file path is shorter, but the problem still isn't solved unfortunately..
I Tried merging two Raster Layers (Digital surface models), but that's the output I get. Does anyone know what the problem is?
QGIS Version: 3.16.4-Hannover
Input parameters:
{ 'DATA_TYPE' : 5, 'EXTRA' : '', 'INPUT' : ['C:/Users/jonas/Documents/UMEP_energyyearroof/0_alle/Energyyearroof_10.tif','C:/Users/jonas/Documents/UMEP_energyyearroof/0_alle/Energyyearroof_11.tif'], 'NODATA_INPUT' : None, 'NODATA_OUTPUT' : -9999, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'PCT' : False, 'SEPARATE' : False }

GDAL command:
python3 -m gdal_merge -a_nodata -9999 -ot Float32 -of GTiff -o C:/Users/jonas/AppData/Local/Temp/processing_gYGbnK/4f8b13ee82e04fe493639dcfc8ab471d/OUTPUT.tif --optfile C:/Users/jonas/AppData/Local/Temp/processing_gYGbnK/1e7ec3ef9602404c998dee96f23658ee/mergeInputFiles.txt
GDAL command output:
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.16/apps/qgis-ltr/./python/plugins\processing\algs\gdal\GdalAlgorithm.py", line 138, in processAlgorithm
GdalUtils.runGdal(commands, feedback)
File "C:/PROGRA~1/QGIS3~1.16/apps/qgis-ltr/./python/plugins\processing\algs\gdal\GdalUtils.py", line 108, in runGdal
for line in proc.stdout:
File "C:\PROGRA~1\QGIS3~1.16\apps\Python37\lib\codecs.py", line 322, in decode
(result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0x81 in position 236: invalid start byte

Execution failed after 0.45 seconds```


texts: How to calculate mean elevation in a DEM between a specific interval?

texts:
I have a DEM with different elevations between 100 m & 3828 m.
I just want to calculate mean elevation higher than 1500 m? Can I do that?


texts: Locate and delete ArcGIS Desktop Lock files

texts:
I wan´t to simply delete the lock files that ArcGIS writes and doesn´t delete after removing data from the contents window. As my ArcGIS Desktop takes around 4 minutes to start, I am getting really annoyed because I wan´t to display in arcgis and at the same time transform the same data with fme. The ESRI Lock files are really getting  on my nerves.

Any tips on how to locate and delete the lock files are most welcome.


texts: Multi Criteria Evaluation implemented via a website

texts:
I'm wondering if there are any websites where spatial multi-criteria evaluation can be performed by a user selecting criteria from their web browser and the result being presented as an online map.

TerrSet does MCE very well, but I would like to have some web application that can replicate what TerrSet does, and present the result as an online map. Ideally, I'd like to provide a range of spatial data sets that a web user can select from, then rank in importance, then perform MCE, which will then return spatial features ranked by suitability according to the ranked multiple criteria that the user selected. Any examples of this online?


texts:
I am trying to export terrain modelling data for 3D printing in stp format but can't find a plugin.

DEMto3D is an excellent plugin that does exactly what I need but only available in Windows.

Is there plugins for Mac that does the same thing?
Is there an alternative data export process?


texts:
I have QGIS 1.8 on Centos 6.5 (I cannot update this system).
I have a problem with opening a shapefile layer with units in seconds (created in ERDAS).
QGIS doesn't display and doesn't save document in the "Degrees, Minutes, Seconds" option when I set it.

What am I doing wrong and how should I work with this layers?


texts: Creating points on street line between two intersecting streets in QGIS

texts:
I have a table which contains fields "Street", "From", "To" and "Number of potholes".

"Street" showing the street line where the potholes are
"From" and "To" are street names which define part of the street where
the potholes are
"Number of potholes" is the number of points which I want to place on
defined part of the line (between "From" and "To").

I want to do this using QGIS.

I have an idea to download street lines from OpenStreetMap and to geocode each street and select street lines by location. I am not sure about steps after these.

Is there any better solution?


texts: OpenLayers Degrees to radians -180 -> +180

texts:
I need to use OpenLayers to rotate polygons directly with the mouse. For this I use ol-ext which works perfectly.

My problem:
I store the angles in degrees in the database and I need the conversion to radian when I load my page.
I store them in degrees for practical reasons in my user editing panel.

This equation of radians to degrees works perfectly:
(e.angle * 180 / Math.PI - 180) % 360 + 180
So I try to find its opposite but without success, here are some examples of correspondences:

//    deg  ->    rad
// -178.33 -> -3.112424297335325
//  178.69 -> -3.164490064610729

//  70.31  ->  1.22
//  71.62  -> -5.033133821377763

//  -1.15  -> -0.02002663689020867
//   1.24  ->  0.02167345873189741


texts:
Using QGIS 3.4 I am currently creating new and also cloning vector layers with postgres datasource in an existing project using PyQGIS. This works OK but I need to publish these layers and enable update, insert and delete. In other words, do the equivalent of going to Project Properties - QGIS Server, and then check checkboxes for "Published", "Update", "Insert" and "Delete" (in QGIS client software).

Is this possible using PyQGIS?


texts: arcpy.ListFeatureClasses returns empty list

texts:
I'm working with SDE connection and want to find if a layer exists on a server. Code sample is:

ws = r'Database Connections\Connection to sde-srv.sde'
# second option to access connection file, same result
# ws = r'c:\Users\CurrentUser\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Connection to sde-srv.sde'

arcpy.env.workspace = ws
fcs = []
for fd in arcpy.ListDatasets():
    for fc in arcpy.ListFeatureClasses(feature_dataset = fd):
        fcs.append(fc)
print fcs  # []

The problem is that everything is OK if I run this in ArcMap console (fcs contains all the FCs from server), but when runnig it in a standalone script outside the ArcMap app SDE seems to be empty.

On top of it running this script with x32 Python crashes with Win error. With x64 it just returns empty list.

Listing any other sources (GDBs and so on) works well in any case.

arcpy.Exists(ws) returns True in both cases.

*UPD. It's not the same situation as described in 'possible duplicate' question (and few others that I've checked before with almost the same title). There is no problem with full/relative path, absence of FC/FD, encodings, priveleges ...

*UPD2 Win10, ArcGIS 10.2.2.3552, both x32 and x64 (Background processing) Python interpreters installed.
RDBMS: Oracle 11g (seems to be 11.2.0)


texts:
I want to download a raster from a WMS server with Python and OWSLib. I have written this code from an example in a manual.

from owslib.wms import WebMapService
wms = WebMapService('http://mapy.geoportal.gov.pl:80/wss/servic/img/guest/ORTO/MapServer/WMSServer', version='1.1.1')
wms.identification.type
wms.identification.title
list(wms.contents)
wms['Raster'].title
wms['Raster'].queryable
wms['Raster'].opaque
wms['Raster'].boundingBoxWGS84
wms['Raster'].crsOptions
wms['Raster'].styles
[op.name for op in wms.operations]
wms.getOperationByName('GetMap').methods
wms.getOperationByName('GetMap').formatOptions
img = wms.getmap(   layers=['Raster'],
                styles=['default'],
                 srs='EPSG:2178',
                 bbox=(708542, 460329, 731172, 493289),
                 size=(32960, 22360),
                 format='image/tiff'
                           )
                           out = open('polandtest.tiff', 'wb')
out.write(img.read())
out.close()

This code returns error:

owslib.util.ServiceException: Parameter 'width' contains unacceptable   value


texts: ArcMap 10.6 Select Layer by Location ERROR 000732: Selecting Features: Dataset does not exist or is not supported

texts:
I'm trying to calculate the accumulative street length in various jurisdictions in a county. To do this I'm using a search cursor in the city layer and then using the current feature to select streets by location and then sum up the total shape length for that selection and print out the city name and associated shape length. I'm getting an error that ESRI states is due to path name, unsupported type, or data type does not work for this tool. I pasted the path right out of catalog, the city layer is a polygon, the streets are polylines. Here is my code so far:

import arcpy

places = r"C:\Users\jsobryan\Desktop\ArcMap\Scrap.gdb\City_Township_Unincorp"

streets = r"C:\Users\jsobryan\Desktop\ArcMap\Scrap.gdb\Streets"

arcpy.MakeFeatureLayer_management(streets,"lyr")

with arcpy.da.SearchCursor(places,["DIST_NAME"]) as cursor:
    for row in cursor:
        arcpy.SelectLayerByLocation_management("lyr","HAVE_THEIR_CENTER_IN",row[0],"","NEW_SELECTION")
        sum = arcpy.Statistics_analysis("lyr",r"C:\Users\jsobryan\Desktop\ArcMap\Scrap.gdb\table_1",[["Shape.STLength()","SUM"]])
        print row[0] + ": " + sum

Here is the error:

ERROR 000732: Selecting Features: Dataset CITY OF ANYWHERE_USA does not exist or is not supported

Thinking maybe I need to do a make feature layer for cities and and use search cursor to select row and then opt to use selected features in select layer by location parms?


texts: How to install MrSID support for GDAL on Ubuntu

texts:
I'm trying to convert some MrSID files to GeoTiff format. They are mostly over 20GB each. It seems the best way to do this is to use gdal_translate with the MrSID drivers. Unfortunately, I can't find anything past about 2013 describing how to install them on Ubuntu. All answers seem to lead back to this question: How to make gdal_translate support MrSID driver however the question is related to converting to MrSID (not FROM MrSID). There are several tutorials from a long time ago but most are 7+ years old and no longer work:

https://trac.osgeo.org/ubuntugis/wiki/TutorialMrSid

I've tried building GDAL from source but am having trouble getting it to recognize all required libraries:

https://trac.osgeo.org/gdal/wiki/MrSID
https://trac.osgeo.org/gdal/wiki/BuildingOnUnix

I also have some concern with this as I don't want to break the package-installed version of GDAL on my system used by QGIS and several other programs.

There are several tutorials using QGIS, however I can't find any new than about 10 years old. I know people have to work with the MrSID format and anything newer involves the same installation of the driver in GDAL. Short of using Windows, is there any solution out for converting MrSID to other formats?

Thanks!


texts: Remove TileStache cache for a single feature

texts:
In a large parcel map, features are edited and need to be re-styled. I want to keep as much of the cache as possible, but invalidate the cache for features that are edited when the features are saved.

In order to do this, I need some way of determining which tiles (at all zoom levels) contain the feature.

Assuming my features are in PostGIS and served up in Leaflet, what are some ways to determine that intersection?


texts: Calculating elevation at point on line using ArcGIS Pro?

texts:
I have two line layers showing stormwater and sanitary sewer pipes.  The attribute tables for both layers show elevations for the upstream and downstream ends of each pipe.  I have a point layer showing where they cross in planimetric view.  Many of the pipes have more than one crossing point.  The attribute table for the crossing points includes IDs for both pipes.

I need to find out which pipe is higher at the crossing points.  I'm using ArcGIS Pro with a Basic license and unfortunately I have minimal coding skills.

Can anyone suggest an approach for this?


texts:
I am trying to delineate flood-prone areas using QGIS.

I have a vector layer of the river and a DEM of the same area.

I have used the GRASS Module to get watershed by using r.watershed and I have also used r.water.outlet to get the basin.

I don't know what to do next to get the flood-prone area. Any suggestions?


texts: Creating new prior distributions in R-INLA

texts:
I'm working on a project where I'd like to compare the performance of BYM latent models as run in three different software packages (winBUGS, CARBayes, and R-INLA) on the same data. Part of this involves trying to create as much comparability between the three models as possible -- such as trying to specify the same prior distributions for parameters in the BYM models.

So far, I've been able to do a decent job of this between winBUGS and CARBayes. However in R-INLA (which I'm new to), I'm running into a problem; I'd like to specify a gamma distribution for a prior, but the R-INLA only has a log-gamma distribution option. I see from here that's it's possible to create one's own prior distributions for use in R-INLA models, but I'm a bit unclear about how to actually do this.

Here's an example of what I have for my R-INLA code at the moment for my BYM model (I'm modeling count data with a Poisson likelihood):

###My model###
form_INLA<-O~1+offset(log(E))+depriv+ f(London_data$idx, model = "bym", graph = adj,
      hyper = list(prec.unstruct = list(prior = "loggamma", param = c(0.001, 0.001)), ###prior specification###
                    prec.spatial = list(prior = "loggamma", param = c(0.1, 0.1)))) ###prior specification###

###Run the model###
    INLA_BYM <- inla(form_INLA, family = "poisson", data = as.data.frame(London_data),
                     control.compute = list(dic = TRUE),
                     control.inla = list(tolerance = 1e-20, h = 1e-08),
                     control.predictor = list(compute = TRUE))

Also, would any of you know if it's possible to set prior distributions in the R-INLA BYM model for parameters other than the random effects terms (i.e. for the intercept term and beta coefficient term)?


texts: Is it possible to negatively buffer a polygon feature so that the output polygon is a certain percentage of the source polygon?

texts:
I am not sure if this is possible, because it seems like a very complex problem.  I have some polygons with attributes that have an attribute indicating a percentage.  I would like to construct a new set of polygons whose area is proportional to their source polygon's area multiplied by its percentage.  So, for example, a polygon with an area of 10 000 and a percentage value of 90% should produce a polygon with an area of 9 000.

Essentially I need to negatively buffer polygon features by a variable distance so that the resultant feature's area is a desired value.

Is this possible in QGIS?


texts: automatically geo- and ortho- rectify raw Landsat data using DEM

texts:
I have a dataset of 20 Landsat scenes (MSS, TM, ETM+) for 20 years (one scene per year). Scenes are in WGS 84, EPSG 32634 coordinate system and I can easily project  = geometrically correct them to S_JTSK, EPSG 5514, using Project raster Tool (Data Management in ArcGIS 10.1).

As my area of interest is in mountaineous region, I would like to orthorectify scenes using DEM (from Shuttle Radar Topography Mission, res ~ 30m) but I don´t know really how to do this. I can´t find my rational polynomial coefficients (RPCs) needed for automatic ortorectification in my data and really I can´t imagine how to do this orthorectification. Any help please?


texts: Generate ESRI grid with xllcenter/yllcenter in header

texts:
I'm trying to generate ARC/INFO ASCII GRID (*.asc) from point feature class using Create TIN - TIN to raster - Raster to ASCII functions. I'm looking for a way to have the output raster file header referring to center of left lower pixel instead of corner. (xllcenter/yllcenter instead of xllcorner/yllcorner).

Is there a way to get this type of header directly from ArcGIS? If not - is it correct if I just rewrite the header from corner to center and change the coordinates to coordinates of center of the left lower pixel (as shown below)?

This is the header from Raster to ASCII tool:

ncols 1000
nrows 1000
xllcorner 1271000
yllcorner 7269000
cellsize 5.0000

This is the header I need:

ncols 1000
nrows 1000
xllcenter 1271002.500
yllcenter 7269002.500
cellsize 5.0000


texts: What is the expected QGIS behavior when adding an ArcGIS Rest vector MapServer layer to QGIS?

texts:
The instructions here: https://hub.qgis.org/wiki/17/Arcgis_rest are straightforward, but I'm getting the same 1000 features regardless of extent (with the 1000 feature outside the viewable map extent in QGIS). My current map extents should only return back a handful of features.  Does QGIS pass in the map envelope as Input Geometryto the ArcGIS REST Service?  Or am I expected to add it?

I am using the method where you take the http://<host>/arcgis/rest/services/<service>/MapServer/<layer_id> and add query?where=objectid+%3D+objectid&outfields=*&f=json to it. And then adding that via Adding URL as a Vector File in QGIS.

Is this expected behavior? And I need to provide additional query parameters to bring down by return features (like input geometry)?

per Mintx suggestions, I checked WireShark and definitely no additional parameters are being sent (is that by design or that a bug).

QGIS version
2.10.1-Pisa
QGIS code revision
d20c5b7


texts: Arcpy Select by attributes IN operator with Tuple

texts:
I have two layers 1. Buffer zones with a string ID
2. the same buffer zones with the String ID plus an extra field call it zone B

*this process is taking place within a loop already, just to give some scope

Objective - extract all of the first layers buffer ID's and use a select by attributes QRY on the second layer to select the first layer zones that are only labeled zone B

What I have done so far:

-extract ID's from first layer store it in a list then tuple

Zonelist = []
for x in gp.SearchCursor(layer1):
    ID = x.ID
    Zonelist.append(ID)
ZoneTuple = tuple(Zonelist)

now I want to construct a QRY for layer2 that selects all the ID's in the tuple and equal zone b

*I am storing the ID's into a tuple to emulate the parenthesis SQL syntax when using the IN operator within Arcmap

layer2QRY = "ID IN {} AND Labelr = 'ZoneB'".format(ZoneTuple)

layer2 = arcpy.MakeFeatureLayer_management(Allzones, "layer2select", layer2QRY)

Update: the QRY works when there are multiple ID's in the tuple. However when there is only one ID in the tuple ex. ('4321',) a comma appears after the one object and it throws an execute error


texts: Why is rasterio getting a segmentation fault trying to read a valid window of data?

texts:
I was trying to run my script and got this output:

orig file has 104424 rows and 161190 cols total
('window (row_start, row_stop), (col_start, col_stop)', ((55150, 72608), (13802, 31386)))
./my_script: line 3: 20326 Segmentation fault      (core dumped) python -O my_main_functions.py $*

I ran through with pdb and found that the segfault was happening on my call to rasterio.open:

with rasterio.open(orig_data_file) as f:
    data = f.read(1, window=window)

The window was valid.  I thought it might be a memory error, so in my debug session I tried doing so with a window of ((55150, 55151), (13802, 13803)) but I still got the segfault.

Why would rasterio segfault trying to read a tiny, valid window?


texts:
Utilizing only ArcMap 10, is there a way to set up a tool that will automatically detect building footprints from aerial imagery and then build a polygon on them?

It doesn't need to be incredibly accurate, but something fast and automated would help. I have no idea where to even start working on such a thing using Python, if its even possible.


texts: Is use of  OpenLayers-2.11 ,GeoExt-1.1 and ext-4.0.7-gpl supported . I have only included ext-all.js in the script src. Should I include others too?


texts: Mosaic overlapping, irregular raster layers produces grid of NaN

texts:
I am trying to use the SAGA - Mosaic operator in QGIS 2.14.0 to mosaic two viewshed raster outputs I have. The viewsheds are irregularly shaped and overlap considerably, thus the need for the Mosaic operator; however, I keep receiving an output grid of NaN values.

Any ideas why this might be occurring?


texts:
I'm loading an external geojson file with the next code:

var sourceTour2016 = new ol.source.Vector({
      url: 'tourdefrance_2016.geojson',
      format: new ol.format.GeoJSON()
  });

  var layerTour2016 = new ol.layer.Vector({
      source: sourceTour2016,
      style: styleFunction
  });
  map.addLayer(layerTour2016); // geojson laag toevoegen aan map

But how is it possible to have the features in an array?
if i'm doing

features = [];
features.push(sourceTour2016.getFeatures());

I receive an empty array...


texts: Defining spatial reference for multiple raster files

texts:
How do I define the spatial reference for a bunch of *.JP2 files?
import arcpy
from arcpy import env

arcpy.env.workspace = "C:/Ortho"
lrs = arcpy.ListRasters("*","JP2")

for rs in lrs:
    print 'Raster = ' + str(rs)
    arcpy.DefineProjection_management(rs,"C:/spatial_reference.prj")

If I then check the properties of the datasets, there is still no spatial reference define. Do I have to write/save it some how?

If I set the spatial reference via the Edit... button (see screenshot) it works fine.

EDIT
Just realized the script works fine! It may has some problems if the files are already open in a different session.


texts: Launch bat. for QGIS 3 python not working

texts:
I'm trying to build a launch.bat for qgis 3 python. I would want to run python via batch file but I get an error 'no module named processing'

Launch.bat:

@ECHO OFF

set OSGEO4W_ROOT=C:\OSGeo4W64

@echo off
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_env.bat"

@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES

set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%

python D:\Testi2\testi.py

pause

Is the problem within those CALL batch files?

I think the problem might be in the python file I'm trying to run:

    from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingException,
                       QgsProcessingOutputNumber,
                       QgsProcessingParameterDistance,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterVectorDestination,
                       QgsProcessingParameterRasterDestination)
import processing

etc.

I also tried adding..

# Prepare the environment
import sys
from qgis.core import QgsApplication
from PyQt5.QtGui import QApplication
app = QApplication([])
QgsApplication.setPrefixPath('C:\OSGeo4W64\bin', True)
QgsApplication.initQgis()

# Prepare processing framework
sys.path.append('C:\OSGeo4W64\apps\qgis\python\plugins')
from processing.core.Processing import Processing
Processing.initialize()

from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingException,
                       QgsProcessingOutputNumber,
                       QgsProcessingParameterDistance,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterVectorDestination,
                       QgsProcessingParameterRasterDestination)
import processing

But it doesn't start qapplication via Pygt5.gtgui


texts: Driver error when converting MapInfo file to Shapefile with ogr2ogr

texts:
I am trying to convert a MapInfo file to an ESRI shapefile using GDAL's ogr2ogr command as described here, but I am getting the following error:

FAILURE:
Unable to open datasource `Iceland2pts.mif' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> a bunch more drivers listed...

Just to be sure I have updated my GDAL install to the latest version (1.10.1) on Ubuntu, so I'm thinking it has to do with the format of the file (due to this), but I just don't know enough about these to know what is going wrong.

Example format of MapInfo file:

VERSION 300
Charset "WindowsLatin1"
COLUMNS 0
DATA
REGION 156
 4
        -22.556955         64.777794
        -22.555510         64.782730
        -22.544556         64.779297
        -22.546009         64.777359
 5
        -22.493431         64.767670
        -22.486217         64.774597
        -22.483826         64.772453
        -22.486826         64.768799
        -22.493431         64.767670

Does the formatting look incorrect, or is something else eschew?

Thanks!

Update: Example of the possible MID header

VERSION 300
Charset "WindowsLatin1"
COORDSYS Earth
PROJECTION 1, 104
COLUMNS 0
DATA
REGION 1
19814
        -23.8252        65.7438
        -23.8289        65.7445
        -23.8355        65.7472
        -23.8421        65.7503
        -23.8477        65.7531


texts: Linking URLs through actions in .url format for QGIS?

texts:
I need to link some DWG files. Normally I'd do that with a simple

cmd /c start ""[% "Link_DWG" %]""

where "Link_DWG" is an attribute containing the URL of the file. My problem is that I need to link to files stored in an ECM, with its own protocol. The ECM's links work when pasted into Internet Explorer (only) or are stored in a .url file, basically like this:

[InternetShortcut]
URL=protocol://yaddayadda

So now I need to somehow create that format from within QGIS. Is this possible?

Also this doesn't work either:

cmd /c start iexplore ""[% "Link_DWG" %]""


texts: How to make water transparent on openstreemap-carto style?

texts:
I have build up manually a tile server on Ubuntu 14.04. LTS, by following switch2osm's guide and have used the openstreetmap-carto style.

Now I intend to add a base layer, underneath the current one, that will display the contour depth lines.

That's why I want to make the water transparent.

How can this be implemented?

I have seen in the scripts/ directory the generate_roads_colours.py. Should I create a new script for water in a similar manner, or one exists?


texts:
I am trying to install the plugin Photo2Shape to my QGIS Pisa with no success.

Does anyone know how to fix this error?

  No module name exifread.

What I need to do is to add geotagged photos to my map that I am working with.

Is this possible with another tool or plugin?


texts: Google Earth Engine: different stretch options

texts:
I'm working on a map view that detects new ice using Sentinel-1 together with Google Earth Engine.

When adding images to Google Earth Engine I can set the min and max values for my layer in numbers but then in the map view I have the possibility to change change these values and also select different stretch options.

Is it possible to choose one of these options instead of entering a max and min view in the code? Right now my code end looks like this:

Map.addLayer(hh2017, {

bands: ["HH"],
max: -5,
min: -19,
opacity: 1,
palette: ['002bff','a4d5ff','979797','080358','000029' ]

}, 'sentinel-1');

I want it to be Stretch: 90% The reason I need to stretch the image each time is that they all look a bit different.


texts: Transform World data from EPSG:4326 to EPSG:3857

texts:
I have tried different open source tools to transform a world data shapefile from EPSG:4326 to EPSG:3857 but I always have an error. I'm not sure, but is it because of Antarctica?

I have tried with GDAL and QGIS and both are unable to do it. Why is it so complicated and how to do it?

The .prj is:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]


texts: Identify raster data for many xy coordinates

texts:
I do not have much experience with GIS. Currently, I have two raster files in ArcMap. I will have to derive data from these rasters for over 5000 xy coordinates. Below I described the approach I am currently using. However, it is very time consuming like this. Anyone who knows a better approach?

use go to xy tool
enter x and y manually
use identify tool and click on the xy
copy values of both rasters to Excel


texts: Do a "3D" extract from a raster stack at specific space and time location in a sf spatial object

texts:
I have a raster stack (or brick) that represent spatial locations as x-y and time as z.  I have a spatial points shapefile (sf object) representing events for which I have the spatial localisation (also the x-y) and the date of the event (also z).  I want to extract the specific raster value associated to the event day.
This is, I think, a relatively frequent problem that I manage to solve with the raster package, the sf package and a bunch of lapply (see my reproducible example below).  However, after looking at my code, I got that feeling that there must be another, more official, way of doing this.  That maybe it was already coded in the raster, sf, stars or another package in R.
Here is an example to show what I want to do and my code to do it:
Preparing the session
library(magrittr)
library(raster)
library(sf)

Making the raster
set.seed(1123)
ras <- lapply(1:7, function(xx) runif(100, xx, xx+1)) %>%
  lapply(matrix, nrow=10, ncol=10) %>%
  lapply(raster, xmn=0, xmx=10, ymn=0, ymx=10) %>%
  do.call("stack", .)

Making the point shapefile
  points <- data.frame(x=runif(10,0,10), y=runif(10,0,10), date=sample(1:7,10, replace = T)) %>%
    st_as_sf(coords = c("x", "y"))

Building a extract function
extract_3d <- function(row, ras){
  row$date %>%
    extract2(ras, .) %>%
    raster::extract(row)
}

Running the function and saving the results
points$res <- points %>%
  split(., 1:nrow(.)) %>%
  sapply(extract_3d, ras=ras)
points
Simple feature collection with 10 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.392242 ymin: 0.2399111 xmax: 9.754538 ymax: 8.89191
epsg (SRID):    NA
proj4string:    NA
   date                   geometry      res
1     1  POINT (8.950553 8.346364) 1.579130
2     5    POINT (6.264524 3.0421) 5.793717
3     1   POINT (9.754538 8.89191) 1.033745
4     4  POINT (2.038216 6.965259) 4.299902
5     2  POINT (2.716971 3.044953) 2.633210
6     5 POINT (1.718511 0.2399111) 5.200416
7     2  POINT (3.215259 4.115244) 2.283648
8     6  POINT (0.392242 2.198866) 6.711109
9     7   POINT (7.32996 5.796152) 7.738792
10    7  POINT (7.956075 4.053148) 7.291688

My code is relatively simple, but it do not scale up to big data sets.  I ran it on ~4 millions date on a tiled Daymet rasters and it's been running for over 3 days on 24 cores and 128 gb of ram.  The idea here is really to do a space-time uneven extract between a raster brick and a spatial point object in R.
Is there a better, more official and efficient way of doing this that my extract_3d function called in a lapply.


texts: Unsnapping two items that have accidentally been snapped in QGIS 3?

texts:
In QGIS 3 how can I "unsnap" two items that have accidentally been snapped?

So far the only way I have found of unsnapping is to delete the vertex and re-add it.


texts:
To calculate the STI one applies the following formula in the raster calculator:
Power(“FlowAccumulationArea_Raster”/22.13,0.6) * Power(Sin(“Slope_Raster”/0.0896),1.3)

Now I have a problem which I don't understand.
When I apply the formula in the raster calculator, the resulting raster shows randomly pixels with no data, while other pixels have been calculated correctly.
My slope and flowaccumulation rasters are consistend and give values for each cell. In my opinion the raster calculator shouldn't pretend that for some of the cells the equation has no solution; which it does by stating "no data".
Anybody has an idea what might be the cause?


texts: Capturing coordinates of map extents in QGIS print composer atlas?

texts:
Here's what I have done so far:

In QGIS 2.10.1 Pisa on Windows 10, I have a set of polygons overlaid onto Google Maps using
the OpenLayers plugin in QGIS.
Using print composer, I have created an atlas which enables me to
show and export each polygon along with the underlying map to a single .jpg or .pdf

Here's what I'm having trouble with:

In print composer, it is possible to see the coordinates of the map extents for each of the atlas items. These are shown in the item properties area for the map I have added to the composition.
I would like to be able to export the coordinates of the map extents, by either including them in a label on the composition itself, or as a standalone dump.

I do not seem to be able to do this using the available expressions.


texts: Is rule-based labelling available for qgis2web plugin?

texts:
Does the qgis2web plugin support multi-attribute, rule-based labeling?

I'm trying to use the qgis2web plugin to export leaflet maps. The labeling on the map I want to export is rule-based. (eg  if x > 0 label with a || x). An export with this style produces a map with no labels.

Here's what it looks like fully functional in QGIS:

If I get rid of the rule-based label, but still use the multi-attibute, multi-line label (eg “X” || ‘\n’ “Y” ), I get a Python export error and no web map is produced.

So in trying to troubleshoot, I changed all the symbology and labels to the simplest terms, and the map successfully exports with labels (but is ugly and doesn't display the data I want). That leads me to believe it is the multi-attribute, rule-based labeling that is the problem.

So in summary, does the qgis2web plugin support these types of labeling schema and I'm just doing it wrong, or is that a feature that's not available?


texts:
I have 2 Layers: Not_Joined and Joined

Not_Joined: Old digitized data that has X attribute in it.
Joined: Spatially joined layer that has the polylines with X attribute in in the table

Each row in "Joined" has a field that says  how many polylines of Not_Joined have been joined to the one element.

Goal: Compare the X attribute in both tables to give accuracy.

Problem: Table length in Joined is 4X bigger than Not_ Joined, so basic for loop will not work here.

Figured out part of it

Query Table from Joined Layer-> Join on Not_Joined based on X.
This seems to work well.


texts:
I'm looking for an algorithm in Python or a combination of arcgis tools to find similar (but not identical) entities between 2 different feature classes.

I attach a picture to the illustration.


texts:
I have downloaded a shapefile of bedrock geology. I need it to be displayed as it should be with the style applied. The supplier will no give me a qml. but i have noticed in the .dbf file the Cyan Magenta yellow and black numbers are stated.

There are hundreds of entries and i dont have the time or patience to go through each one individually.

Is there a way i can apply the colours onto the shapefile (which is polygons and the extracted csv file creates points which have no xy).

QGIS1.8.0 on windows 7.


texts: I have a special MCD43 daily product with 250 m resolution in .hdf format. Whenever I try to open it in the MODIS Re projection Tool, the error i get is "error in module:read hdr file". What should I do? I want to convert these files to .tiff format. Please help!


texts: How to place point information in a .hgt file

texts:
I am a noob when it comes to LiDAR data.  I normally work with las, lasx or las tiles to create golf courses from around the world for a game.  I have recently have had a request to create a golf course in France.  Unfortunately, they do not have LiDAR data in the formats I work with above.  The potential customer steered me to a website that provides .hgt files for France.
I am able to open these files in Global Mapper, but cannot do anything with them because there is no point data in the file.  How can I work with these files so I can eventually convert them to las/laz format?  If I click on classify Ground button in the las tool withing the app I get the message, "You must first select at least one lidar point".  Obviously there are not points in the file.
Is there a way to work with these files for my purpose?


texts: QGIS how to enable spatial indexes on Spatialite table

texts:
I've created a Spatialite layer with all of the waterbodies and watercourses in Saskatchewan for use in QField. Rendering is quite slow on both my older laptop and my field device - a G5 cell phone. When I look at the layer properties I see that the "create spatial index" button is disabled.

I know that Spatialite supports spatial indexes, but it's not clear from what I've read whether the newer virtual spatial index is created automatically or if I need to enable it somehow. If it is truly virtual then I assume it is created on the fly, which would be a bad thing in the case of QField. It's also not clear whether QField will use spatial indexes but that's a question for a different forum.


texts: Error while loading user-defined function in QGIS

texts:
Why does an error occur when I load a user-defined function in QGIS 3.12.2?
The expression:
@qgsfunction(args='auto', group='Custom'):
def getHeight(layId, itemId, feature, parent):
    layout = QgsProject.instance().layoutManager().layoutByName(layId)
    item = layout.itemById(itemId)
    #item.refreshItemSize()
    return item.rect().height()

The error message:
2020-05-16T21:01:46 WARNING Пользовательское выражение getLayoutByName не действительно
Traceback (most recent call last):
File "C:/PROGRA1/QGIS31.12/apps/qgis/./python\qgis\user.py", line 46, in load_user_expressions
import("expressions.{0}".format(name), locals(), globals())
File "C:/PROGRA1/QGIS31.12/apps/qgis/./python\qgis\utils.py", line 744, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "C:/Users/genf/AppData/Roaming/QGIS/QGIS3\profiles\default/python\expressions\getLayoutByName.py", line 1, in
@qgsfunction(args='auto', group='Custom')
NameError: name 'qgsfunction' is not defined

2020-05-16T21:01:46 WARNING Пользовательское выражение getTop не действительно
Traceback (most recent call last):
File "C:/PROGRA1/QGIS31.12/apps/qgis/./python\qgis\user.py", line 46, in load_user_expressions
import("expressions.{0}".format(name), locals(), globals())
File "C:/PROGRA1/QGIS31.12/apps/qgis/./python\qgis\utils.py", line 744, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
File "C:/Users/genf/AppData/Roaming/QGIS/QGIS3\profiles\default/python\expressions\getTop.py", line 1, in
@qgsfunction(args='auto', group='Custom')
NameError: name 'qgsfunction' is not defined


texts:
I have developed an ArcGIS Desktop AddIn for ArcGIS 10.1 using Visual Studio .Net 2010. The AddIn is established, widely used and by all accounts works well. Many folks have also been using it on ArcGIS 10.2 without issue. The one exception is that I call the AggregatePoints geoprocessing tool and this fails on 10.2 with the message:

Could not load the file or assembly 'ESRI.ArcGIS.CartographyTools, Version=10.1.0.0...

As I said, all other operations, including many geoprocessing calls, are working on ArcGIS 10.2 just fine. Before the tool runs, I check that the user is running an ArcInfo (Advanced) license.

Can someone help why this one tool/DLL is an issue on ArcGIS 10.2 please?


texts: How do I get ArcGIS to show script output as it is executing an FME Workspace or Geoprocessing Tool?

texts:
I have a long running (1-2hrs) FME Workspace that I'm calling through a script. How do get the intermediate results to show in ArcGIS while it's executing?

For example, if I call the FME Workspace directly from the Python console or from the GUI, then the results will show as it is executing in the console or the results window, respectfully. However, it won't show any results when called from the script. The best I can do is show all results at the end using the following.

MyWrapperScript.py

import arcpy

# Load required toolboxes
arcpy.ImportToolbox("MyToolbox.tbx","MyToolbox")

# Execute FME Workspace
arcpy.MyToolbox.MyFmeWorkspace()

# Get all messaged produced by the FME Workspace after execution is finished.
arcpy.AddMessage(arcpy.GetMessages())

Even Geoprocessing tools suffer from the same issue when called within a script.

Is it possible to simply reroute output directly to "AddMessage" as the script runs?


texts: How to Change Image's Color By Changing Its Pixels' Values Based On Some Bands Computation in Google Earth Engine?

texts:
I have this requirement to perform pixel-level computation involving spectral bands of sentinel image.
The code in javascript is as follows:
var ratio_1 = (B04 - 0.35) / 0.215;
var ratio_2 = (B04 - B03) / (B04 + B03);

function some_func(band_value) {
  return Math.max(0, Math.min(1, band_value));
}

if(B10 > 0.2 && ratio_1 > 0.7) {
  var a = 0.3 * (ratio_1 - 1);
  return [0.5 * some_func(B08), 0.5 * someFunc(B03), 0.5 * some_func(B02) + a];
}

if(B11 > 0.1 && ratio_1 > 0 && ratio_2 > 0) {
  var b = 5 * Math.sqrt(ratio_1 * ratio_2);
  return [0.9 * some_func(B08) + b, 0.5 * some_func(B03), 0.5 * some_func(B04)];
}

return [1.8 * B08, 1.7 * B03, 2 * B04];

I want to give each pixel a value based on the returned computed bands' array. So this calculated bands-value will change the pixel color. How to go about achieving for an image of sentinel 2?

UPDATE
I tried doing something like this:
    var b2 = image.select('B02');
    var b3 = image.select('B03');
    var b4 = image.select('B04');
    var b8 = image.select('B08');
    var b10 = image.select('B10');
    var b11 = image.select('B11');

    var ratio_1 = b4.subtract(0.35).divide(0.215);
    var ratio_2 = b4.subtract(b3).divide(b4.add(b3));

    function some_func(band_val) {
      var minVal = ee.Number(1).min(band_val);
      return ee.Number(0).max(minVal);
    }

    if (b10.getInfo() > 0.2 && ratio_1.getInfo() > 0.7) {
      var a = ee.Number(0.3).multiply(ratio_1.subtract(ee.Number(1)));
      return ee.Image([ee.Number(0.5).multiply(some_func(b8)), ee.Number(0.5).multiply(some_func(b3)), ee.Number(0.5).multiply(some_func(b2)).add(a)]);
    }

    if (b11.getInfo() > 0.1 && ratio_1.getInfo() > 0 && ratio_2.getInfo() > 0) {
      var q = ratio_1.multiply(ratio_2);
      var w = ee.Number(5).multiply(q.sqrt());
      return ee.Image([ee.Number(0.9).multiply(csome_func(b4)).add(w), ee.Number(0.5).multiply(some_func(b3)), ee.Number(0.5).multiply(some_func(b4))]);
    }

    return ee.Image([ee.Number(1.8).multiply(b8), ee.Number(1.7).multiply(b3), ee.Number(2).multiply(b4)]);

I tried mapping it to image collection. But I get an error saying function argument not initialized
Am I doing wrong? Please help me on this


texts:
MapInfo v10.5 is used across my organsiation by dozens of users as a thick client installation.  A small group of users have started to have an issue whereby the time taken for a new session of MapInfo to open can last up to 30 minutes.

I have logged onto each users’ machine and experienced the same problem as them so this doesn’t seem to be profile related.

When checking the Task Manager during MapInfo startup, the MapInfo process doesn’t seem to use more than 1% of CPU and struggles to use above 26mb memory.  When I do this on my own machine the CPU usage goes up to 13% and the memory usage quickly rises to 175mb and fully loads in less than a minute.

Anyway, just wondering if anyone else has experienced a similar issue with MapInfo.

Thanks,
Robin


texts:
I have the project to develop a mobile application (Android) being able to get offline maps, go on the field, collect datas about the area we're at and then synchronize it with our data server when we're back to the office.

This is for my end of study work (sry, don't know a best translation in english), and i looked at apps like gvSIG Mini, qGIS for Android, qmap, ..
But none of this really fits my need, as it is hardly adaptable ..

Any ideas on api/existing apps/anything ?


texts: Converting multiple geometry KML to Shapefile

texts:
I have a KML file with multiple geometry, points and polygons. I want to convert it to a point Shapefile and a polygon Shapefile.

Using OSGeo4W Shell I tried this :

ogr2ogr -where "OGR_GEOMETRY='Point'" -f "ESRI Shapefile" D:\env\transit_points.shp D:\env\mosquee.kml
ogr2ogr -where "OGR_GEOMETRY='Polygon'" -f "ESRI Shapefile" D:\env\transit_polygon.shp D:\env\mosquee.kml

I get two empty LINE shapefiles with empty tables with nothing but the row names.


texts: Ingesting Landsat8 Collection 1 Level 1 dataset

texts: I am done with indexing, now trying to ingest. I couldn't find any examples for the ingest config for landsat8 collection 1 level 1, does anybody have any samples out there?


texts: Automatically publishing WMS/WFS in Geoserver

texts: Does any one know how to automate publishing of a WMS or WFS Layer in Geoserver given any data source (vector shape-files) using script or programming language (preferably Python)


texts: Is it possible to get GeoJSON geometry from GetFeatureInfo WMS request from ArcGIS Server 10.x

texts:
I'm trying to perform a GetFeatureInfo request to an ArcGIS server 10.x WMS service and would like to get a response in GeoJSON format including the feature geometry.

Here's the request:
    http://mapservices.gov.yk.ca/arcgis/services/GeoYukon/GY_OilGas/MapServer/WMSServer?REQUEST=GetFeatureInfo&EXCEPTIONS=application%2Fvnd.ogc.se_xml&BBOX=-140.742239%2C64.444492%2C-131.499347%2C69.721576&SERVICE=WMS&INFO_FORMAT=application%2Fgeojson&QUERY_LAYERS=1&FEATURE_COUNT=50&Layers=1&WIDTH=578&HEIGHT=330&styles=&srs=EPSG%3A4326&version=1.1.1&x=201&y=207&

and the response contains a null geometry.

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
       "geometry": null,
       "properties": {
       "OBJECTID": "7",
       "SHAPE": "Polygon",
       "SEDIMENT_BASIN_NAME": "Eagle Plain",
       "SEDIMENT_BASIN_ID": "E",
       "BASIN_LOCATION": "onshore",
       "DESCRIPTION": "Northern Yukon Fold Complex",
       "REGION_NAME": "North Yukon",
       "AREA_KILOMETRES": "20803.8",
       "AREA_ACRES": "5140816.4",
       "AREA_HECTARES": "2080375.7",
       "SHAPE.AREA": "20803756917.106998",
       "SHAPE.LEN": "1391344.050657"
     },
     "layerName": "1"
  }
 ]
}

Is it possible to get the geometry in GeoJSON or another format, perhaps by configuring the ArcGIS Server?