texts:
i want to add raster functionality in django using postgis database using django-raster package and this docs.

here my code :

models.py

from django.contrib.gis.db import models
class RasterWithName(models.Model):
    raster = models.RasterField()
    name = models.TextField()

add raster using console :

from django.contrib.gis.gdal import GDALRaster
from myapp.models import RasterWithName
gdal_raster = GDALRaster('C:/Users/username/Desktop/image.tif')
rast = RasterWithName(name='one', raster=gdal_raster)
rast.save()

but i am very confused with Rendering tiles and urls.

docs say :

  /raster/tiles/layer_id/{z}/{x}/{y}.png where the layer_id is the
  primary key of a raster layer. This structure can be used directly in
  online mapping software such as OpenLayers or Leaflet. An example
  request could look like this: /raster/tiles/23/8/536/143.png,
  returning a tile in png format of the layer with ID pk=23 at zoom
  level z=8 and indexes x=536 and y=143

but indexes x,y mean some random x,y coordinate in image ?anyway i using this url http://127.0.0.1:8000/raster/tiles/1/8/536/143.png or http://127.0.0.1:8000/raster/tiles/1/8/20/40.png(origin x,y in my image) and i take a black window in my browser,any idea where i have wrong in my url rendering ?


texts: Adding information to layers in group (QGIS)

texts:
For a QGIS plugin I want to access the layers in a group called 'Drawings'. For this I want a ComboBox with the list of layers. I tried different methods but can't get past the fact that i am refering to a list. For now this is my code:

from qgis.core import QgsLayerTreeGroup

# Fetch the currently loaded layers
layers = QgsProject.instance().layerTreeRoot().children()

# Access Drawings group and add the childeren to a ComboBox
group = layers.children(QgsProject.instance().mapLayersByName('Drawings'))
self.dlg.CB_newlanduselayer.addItems(group)

Anyone a clue how to enhance this code so i can get a list of the layers in 'Drawings'?


texts:
This is not something I do often, but I recall it being fairly easy. Today however I can't make it work:
Python 3.8.5 (default, Jul 28 2020, 12:59:40)
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import rasterio
>>> from matplotlib import pyplot
>>> src = rasterio.open("ProtectedAreas_152160_mask_5x5.gtiff")
>>> from rasterio.plot import show_hist
>>> show_hist(src, bins=50, lw=0.0, stacked=False, alpha=0.3,histtype='stepfilled', title="Histogram")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/rasterio/plot.py", line 250, in show_hist
    arr = source.read(masked=masked)
  File "rasterio/_io.pyx", line 336, in rasterio._io.DatasetReaderBase.read
MemoryError: Unable to allocate array with shape (1, 87425, 89057) and data type float64
>>>
>>> show_hist(src, title="Histogram")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/rasterio/plot.py", line 250, in show_hist
    arr = source.read(masked=masked)
  File "rasterio/_io.pyx", line 336, in rasterio._io.DatasetReaderBase.read
MemoryError: Unable to allocate array with shape (1, 87425, 89057) and data type float64

Is there something missing here? Some type conversion perhaps?


texts:
I have read a number of posts and tutorials, even the QGIS manual on snapping. I set the snapping options under settings>snapping options; selecting the two layers that I would like snapping on, setting the mode to vertex, the tolerance of 0.00002 map units and checked avoid int. and enable topological editing.

I moved the on polygon on layer one, grabbing the approximate vertex to the adjoining approximate vertex on layer two. at a scale of 1:500 it looks as though it snaps correctly, but when I zoom in to a scale of say 1:1 they are not aligned properly. Am I missing something or attempting to be to accurate?

Thanks for the input.


texts: How can I aggregate coordinate points attributes from multiple similar layers into one layer following a specific route in QGIS

texts:
I'm new in QGIS (and GIS in general). I posted a screenshot to help you understand my question.

I have data from 3 different runs measuring air quality in the streets shown in the screenshot. As you can see, there is some overlap between the runs, and they are not exactly the same, as this is impossible to achieve using a phone GPS.

What I want to do is to aggregate the data based on location, in order to end up with one set of data which follows the route I took during sampling. What I'm thinking is somehow to take one point on the route every, say, 10 meters, and aggregate the data points that are the closest to one specific point based on an average or distance weighted average.

The end result I am looking for is an average of the 3 runs in 1 run.

Is there anyway to achieve this? Or a workaround that would provide something similar?


texts:
I'm on IIS 7, running GeoServer I cannot get GeoServer to work on https and looking for any insight on how to get this to work:

Our website and GeoServer are on same domain
https://mymapapp.com
     http://mymapapp.com:8080/geoserver/

Our website is on https, port 443 w/ certificate installed
GeoServer is on http, port 8080

When a customer views our website they get a insecure content error message because we are pulling tiles from GeoServer over http and not https (http://mymapapp.com:8080/geoserver/)

I'm on IIS 7, running GeoServer through Jetty.

Does anyone know what are steps to configure GeoServer to work on https and still be on the same domain?


texts: Mapproxy tiles are partly black instead of transparent

texts:
I use MapProxy for caching WMS/WMTS tiles and it works fine up to a certain zoom level. But for higher zoom levels MapProxy fills the empty parts of the layers in black. This happens for certain geometries, mostly polygons where I only showe their stroke but they have an empty fill.

Here a comparison of the tiles requested to QGIS server and through MapProxy. First the tile requested directly from QGIS server:

And here the tile requested through MapProxy:

It seems to be a problem with transparencies as the white parts also should be transparent. I define transparent true at the sources configuration in mapproxy.yaml:

sources:
  test_wms:
    type: wms
    req:
      url: http://qgis-server.com/service?
      layers: Test WMS
      transparent: true

I additionally tried (without success) to define a global format and use in in cache:

globals:
  image:
    formats:
      imgformat:
        format: image/png
        mode: P
        transparent: true

Another test was to add the image parameter to the layer source configuration:

image:
  transparent_color: '#ffffff'
  transparent_color_tolerance: 0


texts: Option to add 3-band raster image in read only mode in QGIS 3.4

texts:
I tried to open an 3-band raster image(e.g 1.tif) simultaneously in QGIS 3.14 and also
in another software say  Adobe Photoshop  and then edit image in PhotoShop and save it with same name (1.tif)

When tried to save raster image after editing  in Adobe Photoshop, I am getting error message as given below:

  could not save "....\1.tif" because the file is already in use or was left open by another application

Is there any option available in QGIS 3.14 to open an raster image  in read-only mode so that it can be edited & saved in another software.


texts: How to import r.in.gdal from WCS (only region extent)

texts:
Can I import WCS-Data only for the current region by r.in.gdal?

I tried to import raster data from my WCS into GRASS GIS, using r.in.gdal as proposed in
http://grasswiki.osgeo.org/wiki/Global_datasets#OGC_WCS_-_Albedo_example .

Then I added an BBOX definition to the XML file as follows:

<WCS_GDAL>
<ServiceURL>http://localhost/mywcs?</ServiceURL>
<CoverageName>blw.eros_z</CoverageName>
<Timeout>90</Timeout>
<GetCoverageExtra>&BBOX=2596000,1181000,2597500,1182500</GetCoverageExtra>
</WCS_GDAL>

However, r.in.gdal just starts to load the whole dataset in 2MB chunks, ignoring my BBOX argument. r.in.gdal also lacks a -r option, what would be even better.

Altough, when I try to manually download my region extent with a browser, it works well and I get the wanted GeoTiff file for my extent:

http://localhost/mywcs?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&coverage=blw.eros_z&FORMAT=image/tiff&CRS=EPSG:2056&BBOX=2596000,1181000,2597500,1182500&RESX=2&RESY=2

Background: I have large raster files (12 x 1 GB, erosion risk data of CH http://s.geo.admin.ch/bb8ef7c7f) with overlapping areas. Instead of patching everything together with r.patch (did work well for 2 rasters, less for >10...), I want to get only the area of interest per WCS (Mapserver 6.4.0).


texts: Making a population density map from tabulated (.txt) data in R?

texts:
I have California population density data from from the U.S. Census in an Excel file.

The spatial information is provided in ZIP Code Tabulation Areas (ZCTAs).
How do I plot a population density map with that data using R?

This is what I have been doing:

dataset

library(maptools)
library(RColorBrewer)
library(classInt)
library(maps)

## set the working directory.
setwd("../Downloads")

## load the file
pop <- read.table("mre.txt", header=TRUE, row.names=1, sep=",")

#select color palette and the number colors (levels of income) to represent on the map
colors <- brewer.pal(9, "YlOrRd")

#set breaks for the 9 colors
brks<<-classIntervals(pop$POP10, n=9, style="quantile")
brks<- brks$brks

#plot the map
plot(pop, col=colors[findInterval(pop$POP10, brks,all.inside=TRUE)], axes=F)

I also tried this other way:

CaliforniaZTACs<-readShapePoly("[tl_2010_06_zcta510.shp][2]")
plot(CaliforniaZTACs)
CaliforniaZTACs

but could not associate a population density value from my dataset to the location ID of the tl_2010-06-zcta510.shp file.


texts: Adding a WFS with ajax from Geoserver using Openlayers3

texts:
I am trying to load vector with the following code:

 var sourceWFS = new ol.source.Vector({
      loader: function(extent) {
       $.ajax('http://************/geoserver/Viva/ows', {
            type: 'GET',
            data: {
            service: 'WFS',
            version: '1.1.0',
            request: 'GetFeature',
            typename: 'Viva:building',
            srsname: 'EPSG:4326',
            bbox: extent.join(',') + ',EPSG:4326'

          }
        }).done(function(response) {
          sourceWFS.addFeatures(formatWFS.readFeatures(response));
        });
       },

I don`t seem to be able to get any features. I have a suspicion that I am not setting the typename right. The data is points of buildings.

The url in the geoserver preview layer is the following:

http://************/geoserver/Viva/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Viva:building&maxFeatures=50&outputFormat=text.

I have tested loading a GeoJson layer and it works without fail from the same server and data by using the url-method.

The end goal is to add a new point with WFS-t.

I got most of my code from this example:

http://jsfiddle.net/goldrydigital/h7mk0jpt/?utm_source=website&utm_medium=embed&utm_campaign=h7mk0jpt

Update:

I have changed the ol.format.WFS like follow:

var formatWFS = ol.format.WFS({
featureNS: "http://serverURL:8080/geoserver/vivadata",
featureType: "Viva:buildingType"
                                    });

and I get the following error in the console:

Uncaught TypeError: Cannot read property 'writeTransaction' of undefined at transactWFS

This is the block of code for transactWFS:

var transactWFS = function(mode, a) {
      var node;
      switch (mode) {
        case 'insert':
          node = formatWFS.writeTransaction([a], null, null, formatGML);
          break;
        case 'update':
          node = formatWFS.writeTransaction(null, [a], null, formatGML);
          break;
        case 'delete':
          node = formatWFS.writeTransaction(null, null, [a], formatGML);
          break;
      }

and this is my draw function:

var DrawControl = function(opt_options) {
        var options = opt_options || {};
        this.active = options.active || false;
        var draw = new ol.interaction.Draw({
            source: options.source,
            type: 'Point',
            geometryName: 'geom'

        });
        draw.on('drawstart', function(t) {
        buildingLayer.set('selectable',false);
      });

        //Sets building properties when addBuilding-button isclicked
        draw.on('drawend', function(t) {
        var a = t.feature;
        a.setProperties({

            prim_add: document.getElementById('inpPriAddress').value,
            sec_add: document.getElementById('inpSecAddress').value,
            tert_add: document.getElementById('inpTerAddress').value,
            build_type:document.getElementById('inpBuildingType').value,
            build_mat:document.getElementById('inpBuildingMaterial').value,
            build_own: document.getElementById('inpHouseOwner').value,
        })
        transactWFS('insert', a);
        buildingLayer.set('selectable',true);
      });


texts: Resampling method used in QGIS Raster save as?

texts:
I have what seems like a simple question, but I haven't been able to find the answer in the QGIS documentation.

I have a raster that is 250 m resolution, and I want to resample it to 1000 m resolution. This means that for every 16 pixels in my original raster, I will have 1 pixel in my new raster. I would like to retain the mode (most frequent value) of the 16 pixels.

Based on the suggestion from Radim here, I am using this approach:

Load raster into QGIS
Right click on layer name, select 'Save as'
Change the resolution in the dialog box from 250 to 1000:

However, I haven't been able to figure out what technique is used to go from 16 cells (at my original resolution) to 1 cell (at my new resolution). Based on inspecting some specific pixels, it doesn't appear to be either the mean or the mode.

(Note: I would like to use the mode, so if anyone has a better technique, I am all ears)


texts: Use QGIS with script only without opening QGIS first

texts:
I would like to use the QGIS algorithms without opening QGIS.

I would like to use a script to automate processing with a Windows computer from IDLE (python GUI).


texts: Which application generates *.MAG, *.MOB and *.CAR files?

texts:
I have following type of files (I have tried FME too.I think they are DGPS generated file?)

".MAG" , ".MOB" and ".CAR" extension files many more
".SES" ,"ERR", "LOG" several

I have attached files in the link

Or see the Snapshots

Now i want to know:

Which application generates these files?
How to read these files (esp ".MAG", ".MOB" and ".CAR") and
How to generate spatial data from these files (esp ".MAG", ".MOB" and ".CAR")


texts: How can I specify certain quantiles for renderer classification?

texts:
For graduated classification in QGIS 2.0.1, I can choose between several statistical methods to classify my data (standard deviation, quantile, natural breaks...) and for each method it is possible to choose the number of classes.

How can I choose the range limits for each class? In other words, choosing the quantile method, can I set each class border to be, for example, the 20th, 40th, 60th and 80th quantile rather than 25th, 50th, 75th and 95th?

It seems that QGIS calculates automatically the class range limits but it is not possible to know the corresponding quantile.


texts: Getting equidistant regions from a map with qgis

texts:
I'm builing a webapp for which I need to get some GIS work. I have a csv file containing data from a city, split by district, with the following relevant elements:

lat, long, geoJson(multipolygon for each district)

What I'd like to do is to split the map into X equidistant regions. I need to do this in order to "rank" each region according to different properties.

I'm using qgis and I tried finding a solution with no success, I'm a begginer in GIS but I know how to program, just in case there's code involved in some solution.


texts: Exporting a region to Drive with an Earth Engine App?

texts:
I found an example using Google App Engine w/ Earth Engine's backend to draw a region on a map and export that region's information to your Google Drive. If you open the link, the example I am referring to is at the bottom of the page (it is called "export-to-drive").
https://developers.google.com/earth-engine/guides/app_engine_examples
I just wanted to know if this is possible with an Earth Engine App? Drawing a polygon or giving the user the ability to upload their own shapefile would be nice to know.


texts: R levelplot: How to shade some pixels out to show statistical significance

texts:
I am posting this also on stackoverflow because the topic is in the middle between GIS and R programming.

I need to include statistical significance in a map in R using rasteVis, which is based on levelplot.

Using the example from this answer, I can write some text on each of the cells, or on a selection of those, as for instance:

library("lattice")
x <- seq(pi/4, 5*pi, length.out=10)
y <- seq(pi/4, 5*pi, length.out=10)
r <- as.vector(sqrt(outer(x^2, y^2, "+")))
grid <- expand.grid(x=x, y=y)
grid$z <- cos(r^2)*exp(-r/(pi^3))

p <- levelplot(z~x*y, grid,
               panel=function(...) {
                 arg <- list(...)
                 panel.levelplot(...)
                 panel.text(arg$x, arg$y, ifelse(arg$z>=0,round(arg$z,1),""))})
print(p)

Instead of writing something, how can I shade out these pixels (using for instance light grey lines)?
Should I overlay a second raster?


texts: WFS filter encoding with multi word values in QGIS

texts:
I want to create a layer with WFS in QGIS, and filter out some objects with a filter encoding. My problem is that if I try to filter out based on a attribute value with two or more words in it, I get no results. I have tried the same url and FE in an web browser, there i do get results.

Does anyone know how to fiks this?

Here is an example that gives result in web browser, but not in QGIS:

url = "http://wfs.geonorge.no/skwms1/wfs.tilgjengelighet_tettsted?service=WFS&request=GetFeature&version=2.0.0&srsName=urn:ogc:def:crs:EPSG::4258&typeNames=app:TettstedInngangBygg&"

FE = "FILTER=<fes:Filter><fes:PropertyIsEqualTo><fes:ValueReference>app:tilgjengvurderingRulleAuto</fes:ValueReference><fes:Literal>Delvis tilgjengelig</fes:Literal></fes:PropertyIsEqualTo></fes:Filter>"

layer = QgsVectorLayer(url + FE, "test", "ogr")

print(layer.isValid())

The value in this examble is "Delvis tilgjengelig", if it is changed to the singel word "Tilgjengelig", it works in QGIS


texts: Invalid Topology: Cannot Open Cache File when trying to do a dissolve

texts:
I'm trying to perform a dissolve to take some data that is at the Census Block level to the County level. I'm dissolving by the County ID and doing a Statistic Field for various attributes (1 SUM, 2 MEANS, 1 COUNT).

However, it keeps failing with this error.

It's a fairly large data-set, since its census blocks.

I'm working in ArcMap.

Any help?


texts: How do I query based on two tiff files in QGIS?

texts:
I'm using two raster (tiff) files of the Global Human Settlement dataset. One is representing the population in each cell, and the other one is representing the type of settlement (1=rural,2=low density urban cluster, 3=high density urban cluster).

I'm trying to create a table with the total population for each type of settlement. Is there a way to create this table using QGIS?


texts: How to redraw / refresh map on canvas resize?

texts:
I've created a web page that contains a button to set the map full screen and back to its normal form, the problem map doesn't set to full screen until I resize the output window, how to overcome this,

Any guidance and help is appreciated.

Here is the working example: Demo


texts:
I'm trying to create a shapefile layer based on the spatial projection of an already existing raster. (Ultimately, this will be used in gdal.Polygonize to get a shape from my raster data.) However, I am encountering an error when trying to use CreateLayer.

raster_uri- location of the raster file that I want to make into a shapefile. I know this exists.
out_uri- Location I want the new shapefile to print to.
layer_name- String layer name input.
filed_name- String field name input.

This is the code that I was using.

raster = gdal.Open(raster_uri)
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.CreateDataSource(out_uri)

spat_ref = osr.SpatialReference()
proj = raster.GetProjectionRef()
spat_ref.ImportFromWkt(proj)

layer = ds.CreateLayer(layer_name, spat_ref, ogr.wkbPolygon)

It yields this error message:

TypeError: in method 'DataSource_CreateLayer', argument 2 of type 'char const *'

Outputing the type of spat_ref shows up as <class 'osgeo.osr.SpatialReference'> so I'm unsure what char const * it's complaining about.

Interestingly, if I try to do this all from a python shell using the same file, layer creates just fine. Any ideas on what might be the issue?


texts: overlay div overrides another in Openlayers 3.9.0

texts:
I have years to play with Openlayers and I am new to Openlayers 3.

I have the following code

var layer = new ol.layer.Tile({
        source: new ol.source.OSM()
    });

var center = ol.proj.transform([-1.812, 52.443], 'EPSG:4326', 'EPSG:3857');

var overlay = new ol.Overlay({
    element: document.getElementById('overlay'),
    positioning: 'bottom-center'
});

var overlay2 = new ol.Overlay({
    position: center,
    element: document.getElementById('overlay2')
});

var view = new ol.View({
    center: center,
    zoom: 6
});

var map = new ol.Map({
    target: 'map',
    layers: [layer],
    overlays: [overlay2],
    view: view
});

map.on('click', function(event) {
    var coord = event.coordinate;
    var degrees = ol.proj.transform(coord, 'EPSG:3857', 'EPSG:4326')
    var hdms = ol.coordinate.toStringHDMS(degrees);
    var element = overlay.getElement();
    element.innerHTML = hdms;
    overlay.setPosition(coord);
    map.addOverlay(overlay);
});

overlay2 should just stay on the top of the map and do nothing.

overlay creates an info box over the map, where you clicked, showing info about degrees etc.

overlay works, but I cannot see overlay2 and I cannot understand why. Console gives no errors and when I hit ol.Overlay.length , says 1.

Does it get overridden?

And why?


texts: Selecting and zooming in on features from standalone PyQGIS script?

texts:
I want to select and zoom in features in active QGIS project from external application.

This script do all the things I need, except it does not affect currently running project. It selects and makes zoom to selected but not in opened QGIS project. It shows standalone python window instead of project's.

If I remove last command - qgis.exec_() I see only empty console window - no changes in running project's window.

# coding=UTF-8
from qgis.core import *
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsMapLayerRegistry
from qgis.core import QgsMapLayer
from qgis.gui import QgsMapCanvas,QgsLayerTreeMapCanvasBridge
import os.path

qgis = QgsApplication([], True)
qgis.initQgis()

project = QgsProject.instance()
canvas = QgsMapCanvas()
if project.fileName() == "":
    my_file = os.path.isfile(os.path.join('c:\qgis_test','test.qgs'))
    if my_file:
      project.read(QFileInfo(os.path.join('c:\qgis_test','test.qgs')))

bridge = QgsLayerTreeMapCanvasBridge(project.layerTreeRoot(), canvas, qgis)
layers = QgsMapLayerRegistry.instance().mapLayers()
expr = QgsExpression( "\"name\"='{}'".format('Gothenburg'))

selLayer = layers.items()[0][1]

it = selLayer.getFeatures(QgsFeatureRequest(expr))
ids = [i.id() for i in it]
selLayer.setSelectedFeatures(ids)
box = selLayer.boundingBoxOfSelected()

canvas.setExtent(box)
canvas.refresh()
canvas.zoomToFeatureIds(selLayer,ids)
canvas.refresh()
canvas.show()
qgis.exec_()

I think I need to create link between QgsApplication in script and currently running project window but don't know how to do it.


texts: Postgis - st_MakeLine and st_lineLocatePoint

texts:
Lets explain, I have 4 points (1,2,3,4), and I need to visit clients (A, B, C, D).

So, I need to use st_makeline with (1,2,3,4,3)...

Now I have to get % position from A,B,C,D; but in visit order.

If I do a st_lineLocatePoint, A, B and C, returns ok, but D, return first equal position (lesser that C).

How can I solve this? Some ideas?


texts: How to make a copy of a temp raster in QGIS?

texts: I have a temp raster (SAGA raster) file. How can I make a copy of this file (like "save as")?


texts: How to set default language in GeoNetwork 3.0?

texts:
I have GeoNetwork 3.0.4 running and every time I go to the address, it redirects me to the Portuguese version

geonetwork/srv/por/catalog.search#/home

I'd like for all users to load only the site in the English language

geonetwork/srv/eng/catalog.search#/home

Where are those settings controlled? I tried changing the settings in config.xml but that did not work.


texts: How to increase Java heap size on GeoServer on Windows?

texts:
I tried setting wrapper.java.maxmemory=1300 parameter in c:\Program Files (x86)\GeoServer 2.11.1\wrapper\wrapper.conf but I can't set more than 1300MB (geoserver service can't start).
I tried configure Java with -Xmx2048m, for example, command line argument but it ignored.
How to increase max heap size up to 8GB, for example?
Problem is than the GeoServer can't generate big rasters or kmz's on GetMap operation. Maximum about 500MB.


texts: Preventing Append_management from adding target to Table of Contents?

texts:
I am using ArcGIS Desktop 10.2.

I have a script that successfully converts a .csv into an already-existing feature class by first truncating the feature class table and then using the following line to replenish the table with the new data from the first parameter, "source" (an XYEventLayer).

arcpy.Append_management(input, target, "NO_TEST", mapping)

The destination feature class ("target") already exists in the TOC, within two group layers.

But every time I run the script a new duplicate layer appears at the top of the Table of Contents.  I'm 99% sure this is being done by the Append_management tool although ESRI's documentation doesn't mention this behavior.

How can I prevent the Append_management tool from adding its output ("target") to the TOC?


texts: Using getValues in R but unable to obtain corresponding row number and column number

texts:
I currently have a list of PC values for each row and column of a raster. I would like an output that looks like this:

x     y     PC
1    500    -3.45
1    501    -2.8

Here, x and y stand for row number and column number.

I am using the current code in R, but I end up obtaining a matrix only.

r <- pre_pca$map$PC1
r
class       : RasterLayer
dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : in memory
names       : PC1
values      : -9.283828, 3.889567  (min, max)

s <- getValues(r, format="matrix")
head(s)

  X1       X2       X3       X4       X5       X6       X7       X8       X9      X10      X11
 1 NA       NA       NA 3.121756 3.056287 3.030809       NA 2.977990 2.954814 2.869022 2.802258
 2 NA       NA 3.154253 3.116532 3.023211 2.984813 2.973805 2.917697 2.896259 2.795462 2.754197
 3 NA       NA 3.079315 3.038321 2.985467 2.923184 2.909316 2.853848 2.858537 2.767026 2.749681
 4 NA       NA 3.083964 2.983488 2.948048 2.915955 2.879146 2.859355 2.748576 2.760476 2.724777

I have shown only the first 4 rows and x number of columns, but I would not want an output like this. In addition, please help removing NAs as well.


texts: Create custom projection for Shapefile, given specific target coordinates

texts:
I have a somewhat odd situtaion. A client has provided me with a Shapefile that has no projection information, along with an image of the map that it corresponds to. I have run some coordinates through WhatTheProj and it matches no known EPSG code.

So, my next thought is to somehow create a projection that maps coordinates correctly. I know how to create a Proj4 string and test whether it's right, but I don't actually know anything about projections or how to go about this. Are there any tools that would help?


texts: Command-line tool for writing a print-composer to file?

texts:
I'm looking for a command-line tool, to write a previously generated
print-composer to an image-file.
The following resources I already saw:

Related questions (1, 2)
A gist

Nevertheless, all them are not doing, what I'm intending. I'm looking for a tool which looks like

qgisPrintComposer --print-composer 'my-composer1' --output image.png project.qgs

A starting point is here.


texts: How to change values in data slot of a SpatialPolygonsDataFrame in place

texts:
I have a ESRI shapefile with which I do further work in R. It contains polygons of administrative units and their keys. Key value of one unit is incorrect and I need to rewrite it to be able to join my data later on (it is not practical to change the key value in my data).

Changing the key value in the sp object directly gives me an error:

srcShape[srcShape$KOD_ORP == "1000", ]$KOD_ORP <- "3018"

Error in as.vector(data) :
  no method for coercing this S4 class to a vector

I am able to sidestep the error by first taking the entire data slot out of the sp object, changing it as data frame, and pasting it back in.

asdf <- srcShape@data
asdf[asdf$KOD_ORP == "1000", ]$KOD_ORP <- "3018"
srcShape@data <- asdf

I am however not comfortable rewriting the data slot in such a blunt way; it contains hundreds of polygons and I am afraid of messing up the polygon / key relationship.

Is there a safer solution for this? And is the behavior I see a bug, or a feature of Spatial DataFrame objects and exists there for a reason?


texts:
I am using gdalwarp in python to extract and re-project reflectance bands from MODIS level 1B (MOD021km). An example of the command that I am using is the following:

gdalwarp -of GTIFF -tps -t_srs EPSG:4326 -tr 0.01 0.01 -r bilinear -overwrite HDF4_EOS:EOS_SWATH:"C:\tmp\MYD021KM.A2007094.0935.006.2012075055042.hdf":MODIS_SWATH_Type_L1B:EV_500_Aggr1km_RefSB C:\tmp\Output.tif

It converts the hdf subdataset without any problem but the output file seems to have some sort of artifact, specially close to the edges. (See below)

I converted the same data using Modis Conversion Toolkit (MCTK) and it seems to be much better. (Below)

I tried different resampling method and different output projection but still the same results! Is there anyone having this issue with gdalwarp?

P.S. The data can be downloaded from this address:
ftp://ladsweb.nascom.nasa.gov/allData/6/MYD021KM/2007/094/MYD021KM.A2007094.0935.006.2012075055042.hdf


texts:
I have about 600 GeoTIFF map files, each about 140 mb in size that I'd like to convert to KMZ files so people can download them and use them with Garmin GPS units, Google Earth, etc.

Anyone know an easy way to batch process these, preferably converting the tiff image to a much smaller jpeg image as part of the process?

I've stumbled upon gdal_translate which I guess might be useful optimising the images into something smaller using the "-co COMPRESS=JPEG" option.


texts: Putting layers to new place (another coordinate system) using ArcGIS Desktop?

texts:
I have a lot of feature layers  with annotation in geodatabase. I need put all layers to a new place (another coordinate system).

Firstly I export all layers to shapefile and then put it to new place.

But I can't move annotation feature class layers to new place (very long time process). I use editor and move it with control point. It happens fast but saving time is too long.

Does anyone know how to do it easier?


texts:
I stored several layers on GeoServer. In the near future, users should have access to certain layers, e.g. User_1 should only have access to Layer_1, Layer_2, Layer_3, User_2 should only have access to Layer_2, Layer_3, Layer_4 and so on.

How can I realize this?

I already read the user/roles/group documentation but I haven't found a solution for "multiple users and multiple, but overlapping layers".

Or should I create a single workspace for each user? Does that make any sense?


texts: How do I round up the vertical resolution of ASCII files in QGIS?

texts:
I have a large number of .asc DSM  files which I'd like like to lower the vertical resolution for by two decimal places, preferably using a batch operation. I'm hoping this will reduce storage capacity needed and help speed up geoprocessing. Is there any way I can do this either in QGIS or windows command prompt?

To clarify: As an example, I'd like to go from a vertical reading of 1.161 to 1.2 (random raster cell example).


texts: Keep empty GRID Cells in PostgreSQL (after join)

texts:
I want to join a 100*100m Grid  with a point layer. The points contain the number of people living in the grid cell. There are some cells with no polulation, which I want to keep with the Attribute "pop_count = 0". My SQL-Query unfortunately gives back only the cells in which the attribut >0 exists (see picture).

    DROP TABLE IF EXISTS public.routing_join;
    CREATE TABLE public.routing_join as
    SELECT pop_count, grid.wkb_geometry
      FROM population, grid
      WHERE ST_Intersects(population.wkb_geometry, grid.wkb_geometry);

How can I change my query to keep the empty cells?


texts:
I use osm2pgsql to load osm.pbf Data into a PostGIS database.

One objective is to store the osm_timestamp into a column.
The flag "-x" is used in the cmd code.
With the following line in the mapfeatures.style, the osm_timestamp is recognized:

node,way    osm_timestamp   timestamptz(0) linear

This results in a column "osm_timestamp", but all lines are "NULL".

I can't find a solution, because everything but NULL is expected.


texts: Select tolerance or snapping in OpenLayers 3?

texts:
I'm using openlayers 3 in a mobile  application.

My problem  is that I find a big difficulties to  select features using my finger even with a touch pen (I have to be more precise in my click  or touch, thing which is more painful in my user experience).

I looked for something that can help me like the select tolerance or snapping but I found nothing.

Can someone show me the right way to resolve my problem?


texts: Remove duplicate labels between two identical points

texts: I have two spatially identical points that represent two apartments in one building under one address. I only need to label the address, not the apartment number, and the client does not want me to merge the points into one. So, I'm now left with two identical labels for one building. Is there an easy way to remove duplicate labels that exist between two separate features? I just need to show the address once.


texts: How to create a vector feature drop target?

texts:
I would like to create, inside an openLayers map, a drop target for my features.

The user should be able to drag a vector feature, drop on the target so that some further action can be triggered (in the specific it would be moving the feature to a new map object.)

I have tried to do this by adding a OpenLayers.Control.Buttonto the map, my problem now is that I cannot sort out the z-index issue: when dragged on the target, the feature disappears behind the button.

Is there any way i can set up properly the CSS of the button so that it stays above the map tiles but below the features?

I have noticed that the map is on z-index 749, the vector features in a nested div with z-index 330 - I have tried to set the button to z-index 749 too but the feature still moves below.

Ideas?


texts: Network Analyst in Arcgis to generate the network for further analysis what should the attribute table must contain?

texts: I was trying to use the function Network Analyst in Arcgis to generate the network to display the small world phenomenon, so i have to transport the shape file for winbugs to use it for further geostatical analysis, so I was wondering that what should the attribute table must contain in the shapefile to let me make an adjacency matric for further calculate.(my shapefile type is polyline or point)


texts: Unable to connect PostgreSQL 9.5 from QGIS 2.14.1

texts:
I created new database called "shashi" in PgAdmin III. In QGIS to create new connection I entered following details.

In PostgreSQL 9.5
pgAdmin lll
Please enter password for user postgres on server PostgreSQL 9.5: ****

In QGIS
Create a New PostGIS Connection
Name: shashidhar
Service:
Host: localhost
Port: 5432
Database: postgres
SSL mode: disable

Authentication
Username: postgres
Password: **** (same password that I used in pgAdmin lll to connect to the server )

Test connection ERROR:
Connection failed-consult message log for details


texts: Turn on labels for an ArcGIS dynamic MapServer layer

texts:
While I am generally knowledgeable about using the ArcGIS Server REST API to work with dynamic MapServer layers, this is the first time I have tried to turn on labels.

Below is a call for a tile that will display a single green circle.
The value of POITYPE is Visitor Center.
I would like the label ‘Visitor Center’ to appear near the green circle.

Alas, I have been unable to conjure up just the right syntax to make that happen.
What am I missing?

https://mapservices.nps.gov/arcgis/rest/services/NationalDatasets/NPS_Public_POIs/MapServer/export?bbox=-12205302.24789184,4667760.897033789,-12200410.278081588,4665802.198183983&size=1024,410&dpi=96&format=png24&transparent=true&bboxSR=3857&imageSR=3857&dynamicLayers=[{'source':{'type':'mapLayer','mapLayerId':0},'definitionExpression':"upper(POITYPE)='VISITOR+CENTER'",'drawingInfo':{'renderer':{'type':'simple','symbol':{'type':'esriSMS','style':'esriSMSCircle','color':[0,255,0,255],'size':8,'outline':{'color':[0,0,0,255],'width':1.5}}}},'showLabels':true,'labelingInfo':{'labelPlacement':'esriServerPointLabelPlacementAboveRight','labelExpression':[POITYPE],'useCodedValues':false,'symbol':{color:[0,255,255,255],'size':8,'fontFamily':'Arial','fontStyle':'normal','fontWeight':'bold','fontDecoration':'none'}}}]&layers=show:0&f=image


texts: Creating travel time isochrones using ArcGIS Desktop?

texts:
I am trying to build isochrones of 5, 10, and 15 minutes walking time from service area points using ArcMAP 10.4.1. I used Digimap's Integrated Transport Network merged with Digimap's Urban Path as network source data. The first result of running the analysis, i found that some points do not seems to generate correct isochrone, or even no isochrone at all, as seen below:

I then rebuild the network and for this second time using only the Integrated Transport Network data and exclude the Urban Path, resulting:

I also have another different points dataset but when I tried network analysis to those data with the exact same first method (using ITN and Urban Path merged as network data) I did not experience any problem.

Does anyone know why and how to correct this?


texts:
I know this was already discussed in other questions, but I feel I have to ask again.
I have a point layer (locations) and a line layer (electricity grid) which cover the whole world. I have tried using more than one projection (azymuthal, world mercator) but the results weren't good. Is there any possibility to measure point-to-line distances more or less accurately for the whole world? (the accuracy doesn't need to be 100%).


texts: How to Compact geodatabases in multiple directories using ArcPy?

texts:
I am trying to write a script to walk through all folders in my workspace and find and compact all geodatabases in my workspace.  I feel like I'm on the right track but I am stuck.

import acrpy, os
from arcpy import env

arcpy.env.workspace "C:"

for (path, dirs, files) in os.walk(arcpy.env.workspace)
    for file in files
        if ".gdb" in file:
        workspaces = arcpy.ListWorkspaces("*", "FileGDB")
            for workspace in workspaces:
            arcpy.Compact_management(workspace)

I am using ArcGIS 10.2.2 for Desktop.


texts: Posting still running after nearly 4.5 days

texts: I haven't reconciled and posted in about a month. I added a large amount of new data to our database on many different feature classes over this time. When I reconciled I believe there were around 8000 conflicts and I resolved those in the edit versions favor. This all worked fine. The reconcile is done but when I go to post to Default in Arcmap my window freezes. I am currently running the Reconcile Versions tool in ArcCatalog on another work station. It started running on the Oct 14 2016 and is still running now on Oct 18 2016. I have the task manager window open on the work station and it shows the CPU usage going from 0 or 1 up to 30 or 40 and in between in 5 second blips every 20-30 seconds. I believe that the post is still running but I am concerned about how long this is taking.


texts: Deleting sliver polygons after difference operation in QGIS?

texts:
I'm trying to find the differing polygons between two layers in QGIS. After running the difference operation, however, I get results like this:

I am only interested in the solid (blue) polygons, not the weird lines - which I suspect are the result of rounding error, or previous conversion of the polygon layers from different file types - who knows.

One idea I had was to calculate the area of the polygons and simply select those with non-zero (or above any obvious peak in the distribution) values, but irritatingly the line polygons seem to have substantial areas in many cases.

Any ideas?


texts:
I am trying to create a precipitation grid for my meteorologic model.

The problem is that, I do not know what are the things I should have for me to create the grid.

Can anyone tell me what are the requirements and how to run the process?

I am using ArcMap 10.1.


texts:
I'm drawing lines in QGIS and need to specify the width as map units because I want the width to be the same, regardless of scale. The projection I'm using is standard in my field so I don't want to use a different.

My projection map units are meter and I'm working on continental/global scale. The lines should be 2-300km wide.

Everything works fine, except when I put a line width over 100km, QGIS changes the width down to 99999 (m) automatically. Is there any walk around that allows me to specify broader symbols for lines?

QGIS version
2.18.14


texts: I am using the Edit Widget and the Batch Attribute Editor widget to edit our Feature Layer with a Web App Builder.     I would like to be able to modify two attributes automatically when the Edit is Complete with these widgets without the user having to enter the attributes.    I would like to update a date field and a user field.    I know this can be done via the Editor Tracking Fields,  but the layer is edited in ArcMap and I only want to update the attributes when an edit occurs via WAB.     Is there a place in WAB where I can put some custom code to set these attributes when an edit is compete?   Do I need to create a widget that listens for the on edit complete?


texts: How to implement Geography or geometry with spatial index implementation to find entities in the bounding box?

texts:
I have entities saved in SQL table with Lat and long values. I added another column to store Geometry values.

Code for updating Geometry column
UPDATE [dbo].[table]
SET [GeoLocation] = geography::STGeomFromText('POINT(' + CAST([Longitude] AS VARCHAR(20)) + ' ' +
                    CAST([Latitude] AS VARCHAR(20)) + ')', 4326)

Now i need to modify my current procedure which work on Longitude & Latitude columns to newly added Geometry column.
The stored procedure basically takes input as @Minlat, @MaxLat, @MinLong & @MaxLong and returns the entities between these coordinates.These points are then displayed on UI

Current Procedure filter
select Entity, *
from Table
WHERE Latitude BETWEEN @Min_lat AND @Max_Lat
AND Longitude BETWEEN @Min_long AND @Max_Long

I need use the Geometry column in the above where condition so that it returns the same result.

I have created a polygon with the 4 inputs and used STIntersects function however it is not giving me any result. Not sure i should be using Geometry or geography time.
The entities i show on map are static.

Code to create bounding box **
DECLARE @boundingRect Geometry
SET @boundingRect = 'POLYGON((' + @minX + ' ' + @minY + ', ' +
                                  @maxX + ' ' + @minY + ', ' +
                                  @maxX + ' ' + @maxY + ', ' +
                                  @minX + ' ' + @maxY + ', ' +
                                  @minX + ' ' + @minY + '))'
**New Procedure filter
select Entity, *
from Table
WHERE GeoLocation.STWithin(@boundingRect) = 1

If i compare data for the select statments it is different. the STWithin does not select any records from table


texts: Label using VBscript in ArcMap - multiple IF statements?

texts:
I have 4 fields that hold number values.  I only want a label returned if the field <>0.  For instance, if col1<>0, then label = col1 AND if col2 <>0 then label = col2, etc.  So, I want all four columns to show up as labels, but only if they are not equal to 0.
I have tried finding info online and am having issues.  Thought maybe everyone could help.
The following script works, but only for one at a time.  I need to combine them so that this script works with Seq2, Seq3 and Seq4.  I am thinking there is a way to do multiple IF statements?

Function FindLabel ( [Seq1_TractDepthSequence], [Seq1_BegInterval] ,[Seq1_EndInterval], [Seq1_DeckCorpRITot], [Seq1_DeckCorpORRITot] , [Seq1_DeckCorpIntTot] , [Seq1_DeckCorpGWITot]   )
if ( [Seq1_TractDepthSequence] = 1) AND ([Seq1_DeckCorpRITot] <> 0) then
  FindLabel = [Seq1_BegInterval] &" - "& [Seq1_EndInterval] &" : "& [Seq1_DeckCorpRITot] &" RI"
end if
End Function

The following won't work:

Function FindLabel ( [_Seq1_TractDepthSequence], [_Seq1_BegInterval] ,[_Seq1_EndInterval], [_Seq1_DeckCorpRITot], [_Seq1_DeckCorpORRITot] , [_Seq1_DeckCorpIntTot] , [_Seq1_DeckCorpGWITot]   )
if ( [_Seq1_TractDepthSequence] = 1) AND ([_Seq1_DeckCorpRITot] <> 0) then
  FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpRITot] &" RI"
if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpORRITot]  <> 0) then
  FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpORRITot] &" ORRI
if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpIntTot]  <> 0) then
  FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpIntTot] &" CorpInt"
if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpGWITot]  <> 0) then
  FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpGWITot] &" WI"
else
end if
End Function


texts: Is it possible to get a binary geometry from PostGIS?

texts:
PostgreSQL (via C# Npgsql driver) returns geometry column (PostGIS) like a hex string (EWKB as I understand) for this query:

SELECT ST_Transform(geometry, 4326) FROM myGeometries;

Result (for example):

  0101000020E6100000BE6BD0975E595240400E942584A04E40

Is it possible to get it as blob/byte array? Or PostGIS store'em always as text? It's not a good for performance, because you need to convert hex string to a byte array before reading.

PS. I know about ST_AsText, but I don't want to parse WKT, currently I read from EWKB without any problems. Guys, I'm sorry for my primitive English.


texts: Merging polygons which have different coordinate system with map coordinate system?

texts: in order to merge the polygon I go to editor tool and then select merge, however in my case when I want to edit it is giving me this statement" the layers are in a different coordinate system than the current map coordinate system. you can edit data in a different coordinate system than the map" So what can I do so that I can start editing and then activate the merge toolset so that I can merge the polygons. The merge toolset is not active.


texts:
I have some WMS layers with attributes and I want to show some info when I click on some of them.

Currently I'm using:

info = new OpenLayers.Control.WMSGetFeatureInfo({url: 'http://localhost:8080/geoserver/SIMI/wms',
                title: 'Identify features by clicking',
                queryVisible: true,
                eventListeners: {
                    getfeatureinfo: function(event) {
                        map.addPopup(new OpenLayers.Popup.FramedCloud(
                            "chicken",
                            map.getLonLatFromPixel(event.xy),
                            null,
                            event.text,
                            null,
                            true
                        ));
                    }
                }
            });
 map.addControl(info);

 info.activate();

Inside the init() function and works perfectly. When I click something it shows the info.

But I want the popups to be opened as default, when I load the page containing the map. Is this possible? Also, I'd like to be able to open the popup only for a layer (I'm using a layer group and currently using content.ftl to filter out stuff but it's not viable since it still opens a popup with no content).


texts:
I have troubles to get the value of the word file.
I know, the world file includes 6 points.
A, B, C, D, E, F.
I have tried to take a ArcGIS doucument as a reference, however, it does not fit the image to map correctly.
What I have.

Corner coordinates of an image on the screen and on the world map.

Image rotation angle, scaleX, scaleY, skewX, skewY, width, and height of image.

I need the calculation way to find these 6 points correctly.
export const calculateWFContent = (params: any) => {
const {
    mapUpperLeftX, mapUpperLeftY,
    mapUpperRightX, mapUpperRightY,
    mapLowerLeftX, mapLowerLeftY,
    mapLowerRightX, mapLowerRightY,
    imgWidth, imgHeight

     } = params

const A = (mapLowerRightX - mapUpperLeftX) / imgWidth;
const E = (mapLowerRightY - mapUpperLeftY) / imgHeight;
const B =0
const D =0
const C = +mapUpperLeftX + (A * 0.5)
const F = +mapUpperLeftY - (E * 0.5)

return `${A}\n${D}\n${B}\n${E}\n${C}\n${F}`

}
Need to calculate world file content with rotation.
If rotations exist B and D values should not be 0
values:
const params =  {
     imgWidth: '1033',
     imgHeight: '545',
     mapUpperLeftX: '162139.24231791502',
     mapUpperLeftY: '4520936.2860758295',
     mapUpperRightX: '162174.5989677219',
     mapUpperRightY: '4520889.068792953',
     mapLowerRightX: '162149.60850044634',
     mapLowerRightY: '4520870.355746667',
     mapLowerLeftX: '162114.25185063947',
     mapLowerLeftY: '4520917.573029543',
}

Note: Map projection is EPSG:32639


texts: Getting maximum value of a cumulativeCost Image in Earth Engine

texts:
I'm trying to get the max value of a generated cumulativeCost Image in Earth Engine, via the Python API. I thought reduceRegion should work, but it dies.
My relevant code:
malawi = ee.FeatureCollection('USDOS/LSIB/2017').filter(ee.Filter.inList('COUNTRY_NA', ["Malawi"]))
malawi_geometry = malawi.geometry()

nl_raw = ee.ImageCollection('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS')
nl_three_newest = nl_raw.sort('system:time_start', False).limit(3)
nl_stable_lights_only = nl_three_newest.select('stable_lights')
nl_quality_mosaic = nl_stable_lights_only.qualityMosaic('stable_lights')
nl_clipped_to_malawi = nl_quality_mosaic.clip(malawi_geometry)
has_light_at_all = nl_clipped.gte(1)
no_light_at_all = nl_clipped.lt(1)
has_light_at_all_mask = ee.Image().toByte().paint(has_light_at_all.geometry(), 1).updateMask(has_light_at_all)
max_dist = 70 * 1000
dist_to_light_at_all_from_light_mask = no_light_at_all.cumulativeCost(has_light_at_all_mask, max_dist)
dist_to_light_at_all_from_light_mask_max = dist_to_light_at_all_from_light_mask.reduceRegion(
    reducer=ee.Reducer.max(),
    geometry=malawi_geometry,
    scale=30,
    maxPixels=10**9
)
print(dist_to_light_at_all_from_light_mask_max.getInfo())

You can check the whole code with comments and displayed maps in Python here:
https://colab.research.google.com/drive/1z_3_a1UfVwEy1jJEd2Ms1xxw-76HXBPn?usp=sharing
Or in the web editor in Javascript:
https://code.earthengine.google.com/?scriptPath=users%2Fattilastevekopias%2Fdebug%3Anight-lights-debug
The problem is that the final print dies with this message:
HttpError 400 when requesting https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/value:compute?prettyPrint=false&alt=json returned "EEException: Object too large (133984872 bytes)."

What I've tried:

I've checked all the intermediary ImageCollections and Images, put them on a map and they look exactly like they should.
I've tried the exact same reduceRegion function with other images, like the nl_quality_mosaic_max or the nl_clipped_to_malawi_max (the results of the previous steps before the cumulativeCost) and with everything else it works as it should: the print(*.getInfo()) gives back a dictionary like "{'stable_lights': 63}".
I've tried the same code in the web editor (converted it to javascript, see link above) and exactly the same thing happens.
I've tried everything to get any kind of information of the dist_to_light_at_all_from_light_mask_max object like size, type, keys, tried to serialize it and print a slice of that etc etc, but nothing works, I get similar errors, or sometimes just a timeout.
Tried to search for the problem, but google and stack came up with exactly nothing.

So what I'm hoping somebody will be able to help me with:

Why is it not working?
What should I do to make this solution work?
If that's not possible, what's the best way to find out the maximum value in a cumulative_cost image?

The motivation for this that I will generate many similar images in different areas, with different light-intensity cutoffs (or even other types of objects from other ImageCollections) etc, ergo the "maximum distance any point can have in the selected region from a lightsource at least this intense" value will range from hundred to the hundred-thousands meters. And of course I want to display these maps on a colorscale, but to do that I need to know the range it has to cover. And also this exact number is a useful information I need to store for later use.


texts:
The issue:
I am using Windows 7, 64-bit operating system with Python 3.2.  I must continue to use Python 3.2 because of a dependency with another API.  I need some functions in GEOS that are not currently supported by GDAL/OGR 1.8 and there is no port of Shapely for Python 3.2 on 64-bit Windows.

Spatialite has a beta version for 64-bits Windows but I can't get it to work and nobody in the main Stackoverflow list has a solution ( https://stackoverflow.com/questions/8590599/sqlite-load-extension-fail-for-spatialite-in-python ).

I don't want to have to go down the route of using PostGIS to get access occassinoally to one or two functions GDAL/OGR doesn't have....

The Question:
SO:  Can anybody explain how to either get Spatialite working or better still, how I can use CTypes to access the functionality in the GEOS dll (I know how to connect to it but I don't understand how to pass data, say wkt, in and get my result back out - https://stackoverflow.com/questions/8240674/using-geos-via-pythons-ctypes ).

I am VERY sorry for cross-posting but I really do need a solution.  VERY many thanks.


texts: Defining distance between points in meters

texts:
I want to calculate the distance in meters between two points but I am getting
an obscene result 5497980.47338251

Here is the query I made:

SELECT ST_Distance(ST_SeTSRID(ST_MakePoint(6761331.47, 2873883.72), 2154),
                   ST_SeTSRID(ST_MakePoint(2873464.28, 6761335.07), 2154)),
       ST_Distance(ST_MakePoint(6761331.47, 2873883.72),
                   ST_MakePoint(2873464.28, 6761335.07))

The coordinates are in the 2154 SRID.

What should I do to get a good result? I tried with the 4326 SRID and the appropriate coordinates and got a wrong answer as well.


texts: Calculating distance scale factor by latitude for Mercator

texts:
This was really inconvenient topic to google out - so I'm asking for clarification, did I understood things correctly.

I'm talking about measuring short distances in Mercator projection (3857). You can calculate distance in map units (using Pythagoras theorem). In equator, it equals to distance on earth (scale factor = 1); if you move towards poles, distance in map units and in earth are not equal anymore - scale factor starts to increase.

Question is: how to calculate scale factor for given latitude?

Did I understood correctly, formula is factor = 1 / cos(latitude)?


texts: How to create a line between a point and nearest neighbour using python (template included)

texts:
I have 2 shapefiles - one is a series of points, the other is a polyline. I want to draw a straight line between each point to the nearest point on the polyline.

This question has already been asked with regards to ArcGIS. The solution provided by Jeff Berry looks to be exactly what I need:

Creating a line (closest vertex to line)

Unfortunately, I have never worked with python, and have only a basic understanding of coding through Linux.

My questions:

Is anyone willing and able to modify Jeff Berry's solution to fit QGIS?
Or is there an alternative way to create a line between points and polyline, hence avoiding python?*

*NOTE: A useful suggestion is provided here:

QGIS: creating lines from starting to arrival points?

However I've not been able to identify the coordinates of the nearest neighbor - if anyone knows how to do this, I can proceed using the non-python solution.


texts: Where to purchase global boundaries data on city level for commercial use?

texts: Recently I was asked to deal with an app geo locating improve project. It can't be completed without any boundaries metadata on city level. But so far most of my research led me only to things for non-commercial use.


texts: MapserverWCS 2.0.1 GridCoverage, RectifiedGridCoverage, etc

texts:
Actually, I'm working on the implementation of WCS 2.0 service with the latest version of Mapserver (7.0.7).
On GetCapabilities request, the XML response show us a list of layers (wcs:CoverageId) with a CoverageSubtype
for the layer in wcs:Contents, which are available and supported by the server

The first one is,
What is the difference between GridCoverage, RectifiedGridCoverage,
MultiPointCoverage, MultiCurveCoverage, MultiSurfaceCoverage, MultiSolidCoverage, etc ?

Someone can illustrate it with schemas

The second one is,
How Mapserver defines wcs:CoverageSubtype value in XML response from GeoTIFF file ?
In mapfile (*.map), there are no information about grids.

Mapfile example

MAP
    NAME wcs_server
    SIZE 1140 550
    EXTENT -90 -180 90 180
    IMAGECOLOR 255 255 255

    WEB
        IMAGEPATH "/home/tmp/"
        IMAGEURL "/tmp/"
        METADATA
            "wcs_label"                             "Test Mapserver 7"
            "wcs_description"                       "En cours de test"
            "wcs_onlineresource"                    "http//172.17.0.2:9999/cgi-bin/mapserv?"
            "wcs_fees"                              "none"
            "wcs_accessconstraints"                 "none"
            "wcs_keywordlist"                       "wcs test"
            "wcs_address"                           "Rue Truc muche"
            "wcs_city"                              "BLAGNAC"
            "wcs_postcode"                          "31700"
            "wcs_country"                           "FRANCE"
            "wcs_contactelectronicmailaddress"      "test@bla"
            "wcs_contactperson"                     "Quentin SALLES"
            "wcs_contactorganization"               "AKKA-ARTAL"
            "wcs_contactposition"                   "ingenieur"
            "wcs_contactvoicetelephone"             "+33 1 23 45 67 89"
            "wcs_enable_request"                    "*"
            "wms_enable_request"                    "*"
            "wms_title"                             "Test Mapserver 7"
            "wms_srs"                               "EPSG:4326"
            "wms_onlineresource"                    "http://172.17.0.3:9999/cgi-bin/mapserv?"
        END
    END
    PROJECTION
        "init=epsg:4326"
    END

    LAYER
        NAME World
        METADATA
            "wcs_label"                             "Background/World"
            "wcs_rangeset_name"                     "Range 1"
            "wcs_rangeset_label"                    "My label"
        END
        COMPOSITE
                OPACITY 20
        END
        TYPE RASTER
        STATUS ON
        DATA "/home/synext/tmp/Digital_Chart_of_the_World.tiff"
        PROJECTION
            "init=epsg:4326"
        END
    END
END

All layers contain "RectifiedGridCoverage" as value in wcs:CoverageSubtype.
"For an raster layer I can see only "RectifiedGridCoverage" in GetCapabilities and I am not sure if anything else makes sense with MapServer" ---> Where I can found what you said ? In WCS OGC Documentation?
I don't know if it's possible to have others values in wcs:CoverageSubtype ?


texts: While Dragging, created marker with connected polyline from toolbar,previous polyline between two marker deleted using leaflet?

texts:
When dragging created marker with connected polyline from toolbar,previous polyline between two marker automatically deleted  or auto connect from with new draggable marker?

I'm try to dragging marker with polyline,so when i dropping two marker from toolbar and creating polyline between them and then after dragging it is working fine.but while creating another two more marker with polyline and dragging, previous connect marker with polyline,from marker between polyline deleted automatically,and display like that line going to connect/replace with current draggable marker.

so how to solve issue polyline disappear or deleted automatically.

here is code which with i'm trying dragging functionality for marker and polyline after drawing on map from custom toolbar:

var zoomValue = 2;
// center of the map
var center = [42.94033923363183, -108.63281250000001];
// Create the map
var map = L.map('map',{drawControl: false,center,minZoom:1,zoom:zoomValue});

// Set up the OSM layer
L.tileLayer(
'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: 'Data © OpenStreetMap',
    maxZoom: 18
}).addTo(map);

var myFeatureGroup = L.featureGroup().addTo(map).on("click", markerOnClick);

var markers;
var latlngs = Array();
var polyline;

map.on('click', function(e) {
    if(self.isclickedEvent == true){
        var firefoxIcon1 = L.icon({
            iconUrl: self.sourceImage,
            iconSize: [25,41]
           });
        markers = L.marker([e.latlng.lat,e.latlng.lng], {icon: firefoxIcon1,draggable: true}).addTo(myFeatureGroup)
            .bindPopup(customPopup,iDiv)
            .on('dragstart', dragStartHandler)
            .on('drag', dragHandler)
            .on('dragend', dragEndHandler);

        map.fitBounds(myFeatureGroup.getBounds());
        myFeatureGroup.addLayer(markers);
});

function markerOnClick(event){
        if(self.drawPolyline == true){
            latlngs.push(event.latlng);
            if(latlngs.length == 2){
              drawPolyline();
            }
        }
}

//To draw polyline
function drawPolyline(){
    polyline = L.polyline(latlngs, {color: 'blue',weight:6}).bindPopup("PolyLine").addTo(myFeatureGroup);

    // zoom the map to the polyline
    map.fitBounds(polyline.getBounds());
    myFeatureGroup.addLayer(polyline);
    latlngs =Array();
    self.drawPolyline=false;
}

// starting of drag events
var p;
var p2;
var polylineLatlng;
function dragStartHandler (e){
    var latlng = this.getLatLng();
    myFeatureGroup.eachLayer(function (layer) {
    if (layer instanceof L.Polyline)
    {
        p2= layer._latlngs;
        for (var i = 0; i < p2.length; i++) {
            if(latlng.equals(p2[i])){
            p = layer._latlngs;
            polylineLatlng = i;
        }
       }
    }
  });
}
function dragHandler (e){
    var latlng = this.getLatLng();
    myFeatureGroup.eachLayer(function (layer) {
        if (layer instanceof L.Polyline)
        {
            var pipelineLayer = layer;
            p.splice(polylineLatlng,1,latlng);
            pipelineLayer.setLatLngs(p);
            myFeatureGroup.addLayer(pipelineLayer);
        }
    });
}
function dragEndHandler (e) {
    //console.log(e);
    delete this.polylineLatlng;
}


texts: Creating line of varying distance from origin point using Python in ArcGIS Desktop?

texts:
I'm trying to create a line feature from a single point, using a set distance and angle using ArcGIS and Python (ArcPy).

I have a point at these coordinates:  X = 400460.99,  Y = 135836.76

From this point, I want to create a 800 Meter long line at a 15 degree angle from that point.

I do not know what the resulting endpoint will be.

My data are projected in Maryland State Plane South - Meters.


texts: Maximum likelihood classification- changing class colour with every layer

texts:
I created a signature file, where I defined 4 different classes, which I named and coloured specifically.

When I try to use the maximum likelihood classification tool for other layers, the classes defined in the entered signature file have different colours and no names anymore.

It is very important, that the same classes within the different layers can be identified as same classes later on.

Can anybody explain to me, why that is so?

I´m using ArcMap 10.3.1.


texts: Extract from OSM file only nodes with names

texts:
I want to create a Gazetteer for the city I live in. For this purpose I downloaded the relevant .pbf file and converted it to a .o5m file using osmconvert. Now I want to use osmfilter to reduce the file to nodes that are tagged with a name. I am aware that there is a "noname=yes/no" tag for some nodes but that does not cover all nameless nodes.

How can I extract all nodes from the .o5m file where the name tag is given a value?


texts:
I'm running a fairly small script to get height differences. The scripts gets the height in a point based on a DTM. Running the script in the ArcMap python window uses 1 min 30 seconds for 1000 points, while running it in the python window outside Arcmap (C:\Python27\ArcGIS10.5\python.exe) uses 2 min 40 seconds for the same 1000 points. What causes this big difference? There can be many points for each dataset, so I want it to be able to run outside ArcMap so it doesn't run out if memory.

import arcpy
import datetime

arcpy.env.workspace =  r'\\path\Laserdata.gdb'
arcpy.env.scratchWorkspace =  r'\\path\Temp_data.gdb'
# laserline
laser = 'L0380_2m_Vertice_selection'

# 1 x 1 m raster
raster_1m = r'\\path\Connection to sde.sde\DTM1'

def main():
    i = 1
    print "Starttime: " + str(datetime.datetime.now())
    start = datetime.datetime.now()
    fields = arcpy.ListFields(laser)

    if len(arcpy.ListFields(laser,"DTMZ"))>0:
        print "Field exist"
    else:
        print "Does not exist"
        arcpy.AddField_management(laser, "DTMZ", "FLOAT")

    with arcpy.da.UpdateCursor(laser, ['SHAPE@X','SHAPE@Y' ,'DTMZ']) as curs:
        for row in curs:
            X = str(row[0]).replace(".",",")
            Y = str(row[1]).replace(".",",")
            XY = str(X + " " + Y)

            hoydeDTM  = arcpy.GetCellValue_management(raster_1m, XY)
            if str(hoydeDTM) == 'NoData':
                hoydeDTM = -9999
            else:
                hoydeDTM = hoydeDTM

            row[2] = str(hoydeDTM)
            i = i + 1
            if i % 1000 == 0:
                print i
                print datetime.datetime.now() - start
            curs.updateRow(row)

    print "Endtime: "  + str(datetime.datetime.now())
    print  datetime.datetime.now() - start

if __name__ == '__main__':
    main()


texts: Executing ArcMap AddIn function when there is mouse click on map?

texts: I am developing an add-in for ArcMap 10.0 where I need to execute code when there is a click on the map open. How do I do this? I am using VB in Visual Studio 2010. I also need to get the coordinates from the click event.


texts: Uncaught ReferenceError: GEvent is not defined after migrating to google map api v3

texts:
I have managed to change the needed things, but I get the error from the title in the console, and I cannot figure it out what else needs to be done. any ideas? This is my code(if you paste it and run it, it will show):

<!DOCTYPE html>

    <meta name="viewport" content="width=device-width" />
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no">

    <link rel="stylesheet" href="http://openlayers.org/api/theme/default/style.css" type="text/css" />
    <link rel="stylesheet" href="http://openlayers.org/dev/examples/style.css" type="text/css" />
    <link rel="stylesheet" href="http://js.arcgis.com/3.8/js/dojo/dijit/themes/nihilo/nihilo.css">
    <link rel="stylesheet" href="http://js.arcgis.com/3.8/js/dojo/dijit/themes/claro/claro.css">
    <link rel="stylesheet" href="http://js.arcgis.com/3.8/js/esri/css/esri.css">

                <select id="mapServers">
                    <option value="">Please select a MapServer</option>
                    <option

value="http://sampleserver6.arcgisonline.com/arcgis/rest/services/Census/MapServer">Census</option>
                    <option

value="http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/BloomfieldHillsMichigan/Parcels/MapServer">Bloomfield

Parcels</option>
                </select>


texts: How to calculate the scale of map tiles accurately?

texts:
In web maps, there are two coordinate: wgs84 geographic coordinate and web mercator project coordinate. The displayed scale of a world map is the value at equator, and one 256*256 pic can contain the map in web mercator level 0.

So we can calculate the scale like this, 1:(2*pi*6378137)/(256/96*0.0254)=1:591658710.909131. But in the similar Esri provided service is 1:591657527.591555(Esri Ocean Map).

So, what's the affect in map display with the tiny difference? And how Esri calculate his scale?

ps: I was also confused with the same problem in wgs84.


texts: Loading DXF file as vector layer using PyQGIS?

texts:
Can somebody help me with loading a dxf file as a vector layer with PyQGIS?

Is QgsDataSourceURI an option? I've tried it with this code, but it didn't work.

fitxer="/tmp/prova.dxf"
uri = QgsDataSourceURI()
uri.setDatabase(fitxer)
uri.setDataSource('', 'entities', 'OGR_GEOMETRY', "OGR_GEOMETRY='POINT'", 'fid')

layer = QgsVectorLayer(uri.uri(), 'prova', 'ogr')

print 'isValid:', layer.isValid()


texts: Does adding outputs from processing tools in a list save them as layers?

texts:
import processing

layers = qgis.utils.iface.legendInterface().layers()

supply=[]
sub="supply"
for layer in layers:
    if sub in layer.name():
supply.append(processing.runalg("qgis:fixeddistancebuffer",layer,55,10,False,None))

demand=[]
substring="demand"
for layer in layers:
    if substring in layer.name():
        demand.append(processing.runalg("qgis:fixeddistancebuffer",layer,55,10,False,None))

finalints=[]

for s in supply:
    for d in demand:
        intersections_1=[]
        intersections_1.append((processing.runalg("qgis:intersection",s,d,None)))
    finalints.append((intersections_1))
    del(intersections_1)
print(finalints)

Does this save the obtained output files as layers in the list? If so, how can I add the layers saved in the list onto my map? Also, am not able to get intersections of the layers stored in such lists with output as "Unable to run algorithm".


texts: Trying to convert coordinates for use in ArcGIS

texts:
I'm trying to input this data into ArcGIS and plot coordinates like this:

00⁰06.743   031⁰36.970 (coordinates south / coordinates east)

but I've failed so far and I was wondering if someone could tell me how to do this

I've tried to convert these coordinates with http://geoinfo.sdsu.edu/hightech/LM3/dd1.php or GPSvisualizer.com but they only allow you to convert 1 by 1 (I've got about 700 to convert, all saved into Excel spreadsheets)

So I'm hoping for a conversion tool where you can enter a complete spreadsheet at once and get it converted


texts:
I have a mosaic layer that I want to add to my map in ArcMap. When I add this layer manually (i.e the "Add Data" button), a raster image, a boundary, and a footprint are added to my table of contents. However, when I run the following code, only the raster image is added to my TOC. I have no idea why this happens. How do I get my code to also add the footprint and boundary to the map?

    import arcpy
    in_workspace = r"C:\Users\Christian\Documents\ArcGIS\Default.gdb"
    in_mosaicdataset_name = "mosaic_file"
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    path_to_layer = in_workspace+"\\"+in_mosaicdataset_name
    layer = arcpy.mapping.Layer(path_to_layer)
    arcpy.mapping.AddLayer(df,layer,"BOTTOM")


texts: Why can't I buffer polygons in ArcGIS Pro?

texts:
I've created a polygon by merging several polygons from the same layer using the merge option under Edit menu. Now I cannot buffer the resulting polygon. I tried repair geometry it didn't help. I'm getting an error message saying WARNING 000117: Warning empty output generated. How can I solve this? The distance unit I used is meters. I'm using buffer tool. I set side type to "exclude polygon itself", method: planar and no dissolve.

EDIT: I can't buffer unmerged normal polygons either.

CRS EDIT: CRS of the current map: WGS 1984 Aux...
          CRS of the layer to be buffered: WGS 1984 Aux... (set when creating
                                                            the feature class)
          CRS for the output set in the environments tab: WGS 1984 Aux...

NOTE: I'm using ArcGIS Pro 2.2.3.


texts: Searching for Building Footprints or Building Inventory Data for the Balkans

texts: I am searching for free building data and foorprints, preferably GIS shapefiles but data tables are also welcomed, for Albania, Bosnia, and Macedonia. I am working on estimating people per building in a known locations within those countries. Then, applying the values to intrepolate population. Your help and advice is greatly appreciated.


texts: Building on Inside buffer with attributes in ArcGIS?, I'd like to be able to create multiple ring buffers, while retaining (and possibly dissolving on) attribute values. How can this be done?


texts: Find and Replace expression needed in Field Calculator ArcGIS 10

texts: Looking for the Field Calculator code (Python preferred) to replace the last numeral values in a field using another field's values for entire attribute table. Please see screen shot here:


texts: Creating A Legend by Iteration over ee.List() Object Google Earth Engine

texts:
The Goal is to adapt the Legend from the Example for my Data
I just Changed the Following:

In the original the names and their palettes of subclasses are separated in different subclasswise-properties, I guess that is why they need these lines to obtain palette and the names

// Get the list of palette colors and class names from the image.
image.toDictionary().select([BAND_NAME + ".*"]).evaluate(function(result) {
  var palette = result[BAND_NAME + "_class_palette"];
  var names = result[BAND_NAME + "_class_names"];
  loading.style().set('shown', false);

  for (var i = 0; i < names.length; i++) {
    legend.add(makeRow(palette[i], names[i]));
  }
  Map.addLayer(image, {min: 0, max: 17, palette: palette}, 'IGBP classification');
});

the names and palette in the case of CORINE by Copernicus are easier obtained

var names = ee.List(lc2012.get('landcover_class_names'))
var landcoverPalette = ee.List(lc2012.get('landcover_class_palette'))

extracting Values by indexing over theese ee.List() Objects by e.g. names[0] won´t work. Instead one has to use names.get(0)

print()names.get(0) retrieves a Computed Object that looks like a string --> 'Artificial surfaces; urban fabric; continuous urban fabric'

but used inside a function e.g. legend.add(makeRow('E6004D', names.get(0))), it doesnt get the string, but instead delivers this

Code

I do not understand how to access my List and retrieve the wanted element

for (var i = 0; i < names.length; i++) {
    legend.add(makeRow(landcoverPalette[i], names[i]));
}

The legend won´t load fully


texts: Solving QGIS inconsistency between file names

texts:
I'll describe my particular case but I'm sure it's a more general problem. I'm using the r.horizon plugin (GRASS 7.4) in QGIS 3.8.3, which should produce an output raster file for every direction I've specified (from angle 0=East to 180=West, for example). The problem is I don't get all the outputs I should, only the ones for angles >100.
Apparently the software begins the calculations and creates all the files, naming the outputs:

outputda04053511ae49fa9d98db11bc54f4ac_000
outputda04053511ae49fa9d98db11bc54f4ac_090
outputda04053511ae49fa9d98db11bc54f4ac_180

But later when the software is creating layers based on these files, it does a search for:

outputda04053511ae49fa9d98db11bc54f4ac_0
outputda04053511ae49fa9d98db11bc54f4ac_90
outputda04053511ae49fa9d98db11bc54f4ac_180

How do I solve this inconsistency between the name the software gave the files, and the name it uses to retrieve them and create the layers?


texts: Select_analysis in ArcPy gives ERROR 000358?

texts:
Trying to learn how to use Select_analysis but the queries cause the

  ERROR 000358: Invalid expression.

Piece of code:

arcpy.Select_analysis(in_features='test1',
                out_feature_class=os.path.join(r'C:\Users\users\Desktop\05\Project1.gdb',
                 'WACDE'),where_clause="Select test1.KA FROM test1 WHERE test1.NUM = '12'")

What is the correct form?


texts: Show Cypriot names for Cyprus on OpenStreetMap map

texts:
How can I show only cypriot names (instead of Turkish) in OpenLayers map for the northern Cyprus side?

By default OpenStreetMap shows both Turkish and Cypriot.

https://www.openstreetmap.org/#map=13/35.1278/33.6362&layers=C

In OpenLayers it only shows Turkish names.


texts:
I'm currently creating a new website based on geo-localised informations, which is quite new to me.
I need to display my points of interest according to where the map is, and where the user make it move.
I know I can use the Haversine formula to find POIs around the map's center for the first load.

Now, when the user moves the map to the side, how should I proceed? I though I could just trigger an AJAX load to fetch the POI around the new center, but this would not take into account the already loaded POI, so I'd need to clear the map before adding the new corresponding POI.

I think it'd be better if I can only load the missing POI, as it would probably be more efficient on zoom out.

What is considered best practice in GIS to load data on the fly? What is more memory-efficient?

Also, how can I load POI for the area covered by the screen, but load 20 by 20 to show progressively that informations are being fetched?

I'm using leaflet to display data on map, but I don't know if the number of information displayed as markers take more memory when they are grouped in a cluster.


texts: Generating spatially correlated landscapes (2 spatial layers)?

texts:
I want to generate artificial landscapes that vary in "tree cover" abundance and distribution.

This should be relatively easy to do by varying Hurst exponents (level of aggregation) and adjusting the number of pixels representing tree cover in a matrix (in Matlab for example).  The issue is that I also have to include a layer representing enhanced vegetation index. Since EVI is correlated with tree cover in real landscapes, my artificial landscapes should also represent this level of correlation.

I think I can use the Band Collection Statistics tool to determine the correlation between tree cover and EVI of a "real" landscape", but is it then possible to "apply" this correlation value to generate a spatially correlated landscape with a layer representing tree cover and another representing EVI?


texts: QGIS Uploading images or files for features

texts:
I am using a QGIS 2.12 project and I would like to know if it is possible to upload files and photos for each element in a layer.

The files should be stored for example in C:\myproject\files

I can't think of any feature in QGIS right now that could achieve this..

The idea is to use it in a microsoft Surface tablet, so that users can take a picture of a feature with the camera (of a streetlight, a bench..) and assign it to the feature.

Any suggestion?


texts: Why are the results of ee.Reducer.fixedHistogram() fractional (Earth Engine)

texts:
I am calculating a histogram for an area of interest using ee.Reducer.fixedHistogram(). The pixel count per bin is often a fraction, but not always. I was expecting integer results - the count of whole pixels per bin. Why are the results fractions of a pixel?

Here is my code - inspect the histogram table and hover over bars in the histogram chart to see the results:

// Define an image.
var img = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_038029_20180810')
  .select('B6');

// Define an area of interest.
var aoi = ee.Geometry.Polygon(
  [[[-110.5241853160436, 44.82231727808245],
    [-110.5241853160436, 44.474355785401221],
    [-110.06001295276235, 44.47355785401221],
    [-110.06001295276235, 44.82231727808245]]], null, false);

// Show image and AOI on the map.
Map.centerObject(img, 8);
Map.addLayer(img, {min: 100, max: 3500});
Map.addLayer(aoi);

// Calculate a histogram for the pixels in the API.
var swir1HistDict = img.reduceRegion({
  reducer: ee.Reducer.fixedHistogram(0, 5000, 250),
  geometry: aoi,
  scale: 30,
  maxPixels: 1e13
});

// A dictionary was returned, get the ee.Array results.
var swir1HistArray = ee.Array(swir1HistDict.get('B6'));

// Convert the array to a list so that it can be mapped over.
var swir1HistList = swir1HistArray.toList();
print(swir1HistList);

// Convert the array into a feature list.
var swir1HistFtList = swir1HistList.map(function(li) {
  li = ee.List(li);
  return ee.Feature(aoi, {'Value': li.get(0), 'Frequency': li.get(1)});
});

// Convert the feature list into a feature collection.
var swir1HistFtCol = ee.FeatureCollection(swir1HistFtList);

// Chart the histogram.
var chart = ui.Chart.feature.byFeature({
  features: swir1HistFtCol,
  xProperty: 'Value',
  yProperties: 'Frequency'})
  .setChartType('ColumnChart');
print(chart);

Code Editor script


texts: Reverse direction of "Add Geometry Attributes" tool?

texts:
The ArcGIS Desktop tools Add XY Coordinates and Add Geometry Attributes read from [Shape] column and write the shape XY properties to new X and Y columns. I need to go the other direction, from X and Y columns to the [Shape] XY properties. It is a point feature class already, not a bare table.

XY Table to Point isn't applicable because creates a new feature class.


texts:
I want to check Force layer to render as raster using pyqgis:

I have tried:

layer = iface.activeLayer()
layer.renderer.setForceRasterRender #AttributeError: 'bool' object has no attribute 'setForceRasterRender'

and:

renderer = QgsFeatureRenderer()
renderer.setForceRasterRender #AttributeError: 'bool' object has no attribute 'setForceRasterRender'
layer.setRenderer(renderer)

Both generate AttributeErrors. What is the correct way to set it?


texts: Getting Zonal Statistics of Overlapping Circles?

texts:
I have a Buildings (polygon) shp file and a DSM. what I wanted to do was to get the pixel values for all the vertices of those buildings, and find another pixel value just around those vertices (say 2m, but exclude values coinciding with the building) that will give the maximum difference when subtracted to the value of the vertex. Finally, retain only the vertex that gave the maximum difference for each building.

What I did is I used the 'Feature to Point Vertices' in ArcMap to convert my building vertices into a point shp file. Then I made a 2m buffer around these points and I used the Erase tool using the Buildings as the erase feature. I thought I'd get the value for the vertices using 'Extract Values to Point' tool and then use Zonal Stat as Table to get the minimum values on the buffer circles. This is where my problem starts, ArcMap won't perform the Zonal Stat as Table because I have overlapping circles and performing the process on each circle is time consuming.

Does anyone have a suggestion to do it in a more efficient manner?


texts: How to use the HTML  tag together with OpenLayers 3 raster layers?

texts:
In order to include interaction to my OpenLayers 3 raster layer, I would like to use something similar to the HTML <area> tag (http://www.w3schools.com/tags/tag_area.asp). The aim is to click on certain regions on a raster layer (without projection) to trigger an event.

I haven't found anything on the web and would like to hear your experiences, comments etc.

One idea so far is to use an <img> with a transparent gif image (1x1 px) and then extended over the map extent.

<img src="/img/transparent.gif" width="1024" height="1024" usemap="#mapImg">
<map id="mapImg">
    <area id='area1' shape='poly' coords='...'>
    <area id='area2' shape='poly' coords='...'>
    <area id='area3' shape='poly' coords='...'>
</map>


texts: Connect postcode to administrative boundry level

texts:
In OSM I would like depict a certain area ("working area"). Using the following I can specify a certain area around a certain city:

[out:json][timeout:800];
{{geocodeArea:Netherlands}}->.searchArea;
{{radius=40000}}
(
  node["admin_level"="8"](area.searchArea)(around:{{radius}}, {{geocodeCoords:Oirschot}});
  relation["admin_level"="8"](area.searchArea)(around:{{radius}}, {{geocodeCoords:Oirschot}});
);
out body;
>;
out skel qt;

However the radius does not deliver the required area. An alternative is to use the postal code, which might look like:

/*
[out:json][timeout:25];
// fetch area “Eindhoven” to search in
{{geocodeArea:Eindhoven}}->.searchArea;
(
  node["addr:postcode"~"^56"](area.searchArea);
  way["addr:postcode"~"^56"](area.searchArea);
  relation["addr:postcode"~"^56"](area.searchArea);
);
out body;
>;
out skel qt;

I am struggling in finding the administrative border of the latter result, such that I can improve the first code snippet to specify the cities based on a list of postalcodes.

Which direction I should look for?


texts: Rule based styling with select-by-intersection syntax

texts:
I want to do rule based styling for a selection of gridded polygons dependent on their position relating to another layer. I want to show the polygon that a point feature based in an atlas falls within which I have managed with
intersects($geometry,@atlas_geometry)

The tougher bit is that I also want to show the 8 grid boxes that surround the central grid box in another style. Unfortunately, as the grid is rectangular, not square, I cannot just use a buffer of the atlas feature as when the feature is close to the edge it will show features to the east/west before north/south. I have tried using nested intersections such as
intersects($geometry,intersects($geometry,@atlas_geometry))

but no luck so far. I want to do this dynamically so it shows just the grid surrounding the atlas feature and without having to create any more new layers.


texts: Create polygon symbology from attribute?

texts:
I have shapefile with polygons and lines. There are attributes for: layer, lineweight, colour: Red, colour: Green, colour: Blue. (colour values are from 0-255).

Is there a way to read values from attributes, and assign colour and lineweight to polygons and/or lines? After that I would like to use the save layer option to save symbology.

I'm using ArcMap 10.3.


texts: How to add new water bounds to OpenStreetMap?

texts:
Penikkajarvi lake is present on Google Maps but is absent / plain field on OpenStreetMap.

Is the any way (prefer not manual) to add new water body to OpenStreetMap?


texts: How is the Arc coverage "regions" data format represented in shapefiles and geodatabases?

texts:
Prior to shapefiles and geodatabases, the default ESRI vector data storage format was the coverage.  Later iterations of the coverage included the regions concept, whereby (for example) multiple features represented a single region.  Hawaii was a classic example, where multiple polygons (islands) were associated with one region (the cluster of islands comprising the State of Hawaii).  Specific Arc tools were developed to create, manage and analyze regions-based data.  They still exist at the 10.5 level, for example the "Polygon Coverage to Region" tool.

But such tools are specific to the coverage data format, and I don't see the regions concept described in either the shapefile or geodatabase literature.  I'm curious how regions are implemented in those more modern data formats.  Obviously I'm missing something because Hawaii still exists, even with the slow demise of the coverage!

My specific question is, how is the regions concept - as opposed to individual features - stored, managed and analyzed in shapefiles and geodatabases?


texts:
Here is the situation: I've installed QGIS 3.10 with the downloadable .exe file and accepting all as predefined. The Python version used with it is 3.7.0 like here:

import sys
sys.version
'3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]'

The thing now is that I'm working on a script that need GeoPandas and I'm not able to install or import it:

import geopandas
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python\qgis\utils.py", line 737, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'geopandas'

I tried to install anaconda and within it I'm able to install pip and GeoPandas without any issues but it comes with Python 3.7.4. I have no deep understanding about how several Python versions work together within the same PC.

So, my issue here is easy: I need to be able to import GeoPandas in the Python console regarding that I use QGIS 3.10 and it works with Python 3.7.0.

I use Windows 10, by the way.


texts: Converting polyline to presence absence raster in ArcPy?

texts:
I have a polyline that I am converting to a raster in arcpy. This is occuring within a loop and at each iteration I change the processing environment to a different area of interest. In this case it is a specific polygon (like below)

whereClause = '"FID" = ' + str(poly.FID) # Create a clause to select only the current record
arcpy.MakeFeatureLayer_management(shapefile, "currentPolygon", whereClause)
arcpy.env.extent = "currentPolygon"

At each iteration I would like to create a raster where cells containing the polyline are given a value of 1 and cells without polyline are given a 0 within the bounds of the env.extent.

What is the key to creating the 0 values?

Does it come from the arcpy.env.nodata parameter?

I can make the conversion in my loop but cannot quite get the binary format I am envisioning.

arcpy.FeatureToRaster_conversion("roadsLayer", "VALUE", "roadsGrid.tif",)


texts: Create a binary output raster using an algorithm of reflectance indices

texts:
I'm attempting to create a binary output raster in ArcMap based on an algorithm that uses three different reflectance index rasters I've already created from Landsat 5 data (mNDWI, NDVI, and EVI)the identify surface water pixels. The algorithm I need to apply is:

(mNDWI-NDVI>0 OR mNDWI-EVI>0) AND (EVI<0.1)

I need every pixel that the expression is true for to = 1, and every other pixel to = 0.

I tried stacking the three index rasters using the Composite Bands tool and then using the Band Arithmetic raster function, but the algorithm was invalid.

Is there a way to accomplish this in ArcMap?


texts: shapely geometry - how to preserve refrence to original feature

texts:
What is preferred method to save reference from geometry to original feature when working with shapely geometries in python?
Ie. After creating STRtree and getting some result from it, there is no direct way to access feature f from which this geometry was created.
from shapely.geometry import shape
from shapely.strtree import STRtree
...
tree = STRtree([shape(f['geometry']) for f in features])

I guess there is something short easy and cheap as this is often used but can't find it or guess what it is? Dict comes first to mind but geometries are unhashable so?


texts: Removing Trailing Zeros in Geoserver SLD

texts:
I have maps from a PostGIS datastore in Geoserver. The map is styled using SLD created dynamically from my web apps. The problem is that sometimes, the map has float-type attribute data.

The label actually renders correctly using <TextSymbolizer>, but with a trailing zeros. e.g 0.850000000000. What I want is that it should render 0.85 instead of 0.850000000000.

Is there any way I can do that using the SLD?

I'm using Geoserver 2.6.1


texts:
There's a point layer which contains several points without order, I want to create a line to connect them.

I've read similar questions and I can work out with solutions using "point to path" or "point2one" plugins. However, I don't know those points' order unless I see them in the map. So I have to first see the point position in the map, then search the point and edit the order field in the point layer, finally connect them.

Is there any way I can make it easier, imaging just click the two points in the map, it connect and create line immediately?


texts:
I am running into an error running r.slope.aspect I have never seen before.  When I run it, I get the following output:

r.in.gdal complete. Raster map <dem> created.
WARNING: map [slope] - unable to write row 1269
WARNING: map [aspect] - unable to write row 1270
100%
WARNING: error writing null row 1283
WARNING: error writing null row 1284
WARNING: error writing null row 1285

The "error writing null row" continues for a very long time.  I have been running this script many times and only recently am I starting to get this strange error.

echo "Running r.sun for day $DAY"
#Create new projection info
g.proj -c georef=$DEM
#Import Dem
g.mremove -f "*"
r.in.gdal input=$DEM output=dem
#Set Region
g.region -s rast=dem
#Calculate Slope and Aspect
r.slope.aspect elevation=dem slope=slope aspect=aspect
#Create flat map
r.mapcalc "zeros=if(dem>0,0,null())"
r.sun elevin=dem aspin=zeros slopein=zeros day=$DAY step=$STEPSIZE dist=$INTERVAL glob_rad=flat_total_sun
#Using dem and calculated slope and aspect, generate a solar insulation model
r.sun elevin=dem aspin=aspect slopein=slope day=$DAY step=$STEPSIZE dist=$INTERVAL insol_time=hours_sun glob_rad=total_sun
#Output files
r.out.gdal -c input=total_sun output=./global/daily/total_sun_day_${DAY}.tif
r.out.gdal -c input=flat_total_sun output=./global/daily    /flat_total_sun_day_${DAY}.tif
r.out.gdal -c input=hours_sun output=./insol/daily/hours_sun_day_${DAY}.tif
r.out.gdal -c input=slope output=./slope.tif
r.out.gdal -c input=aspect output=./aspect.tif


texts: How to get feature data by passing lat/long value to FeatureServer in ArcGIS Online?

texts: I have created a feature layer and deployed it as a service on ArcGIS Online? Is there any way to pass lat/long values as query parameters and get respective feature layer details?


texts:
The ArcGIS SpatialJoin online help page has a long list of match_options when joining data spatially, including HAVE_THEIR_CENTER_IN.

However running the arcpy.SpatialJoin_analysis tool in Python with this option results in an error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000800: The value is not a member of INTERSECT | CONTAINS | WITHIN | CLOSEST.
Failed to execute (SpatialJoin).

These options are also missing from ArcToolBox:

The licensing on the help page suggests this tool is available for all versions of ArcGIS:

  ArcView: Yes ArcEditor: Yes ArcInfo: Yes

So is this an error in the docs, licensing on keywords, or a misunderstanding on my part?


texts: Looping Sub-folders to merge shapefiles export by sub-folder name to file geodatabase

texts:
trying to loop through folders, merge shapes in folders, naming my folder name, export to feature class in file geodatabase.

**merge loop**

import arcpy
from arcpy import env
import os

**Workspace**

sourceFolder = r'C:\00_TestingScript\Test01'
outputFolder ='rC:\00_TestingScript\Test01_Output'

env.workspace = sourceFolder

folderList = os.listdir(sourceFolder)

try:

    for cntyMems in folderList:
        print(cntyMems)
        folderList = sourceFolder +'/'
        print(folderList)
        arcpy.env.workspace = folderList
        memShapeOutput = outputFolder + "\\MemShap_" + cntyMems

        FCList = arcpy.ListFeatureClasses()
        for f in FCList:
            print(f)
            #########################################################
            #get stuck here
            #arcpy.Merge_management(f,memShapeOutput)
        print("I am working amazing")
except:
    print("sorry")


texts:
I'm making a small thematic map of the Washington DC area using d3.js.

Everything I've tried "tilts" the square of DC.

Is there a d3 projection that will put DC at the right angle?

https://github.com/d3/d3-geo-projection/


texts:
I want to change my text backcolor.

        new ol.style.Style({
        text: new ol.style.Text({
            text: marker_name,
            fill: new ol.style.Fill({
                color: '#fff'
            })
        })
     })

"fill" property change text color.
I search from OpenLayers v4.0.1 ol.style.text but I did not find propery for this.

You think draw a shape under the text like this

 new ol.style.Style({
        image: new ol.style.RegularShape({
            fill: fill,
            //stroke: stroke,
            // offsetY: -25,
            points: 4,
            radius1: 30,
            radius2:10,
            angle: Math.PI / 4
        }),
        text: new ol.style.Text({
            text: marker_name,
            fill: new ol.style.Fill({
                color: '#fff'
            })
        })

But OpenLayers v4 can't draw rectange in style. For this issue you can look that.