texts: Local Nominatim server returns incorrect place with GPS lookup

texts:
I am trying to run a local Nominatim server with a smaller dataset that only includes cities, towns and villages similar to this Stackoverflow question. My goal is to look up which GPS coordinates correlate to what city, I have tried creating a smaller dataset with two different commands:
./osmosis --read-pbf planet-latest.osm.pbf --tf accept-nodes place=city,town,village --tf reject-relations --tf reject-ways --lp --write-pbf planet-excerpt.osm.pbf

and also
./osmosis --read-pbf planet-latest.osm.pbf --tf accept-nodes place=city,town,village --tf accept-relations --tf accept-ways --lp --write-pbf planet-excerpt2.osm.pbf

The problem is that the reverse lookup on my local Nominatim server returns a different place than the official one, I have included some examples of the returned json outputs.
Place 1 local (Chillum) vs Place 1 official (Washington)
Place 2 local (Getafe) vs Place 2 official (Madrid)
Place 3 local (Gatineau) vs Place 3 official (Ottawa)
There seems to be an issue with the correct detection of the bounding boxes, how can I improve this? I'd rather not import the whole planet and if I have to go that route can I drop information afterwards without impacting the reverse lookup?


texts: How to add Domain using ArcGIS Diagramer?

texts: Is it possible to add domain to XML Schema of GDB using ArcGIS Diagramer ?


texts: Where can I find country borders data in plain text format? (XML, JSON, CSV, etc)

texts:
I'm trying to obtain some data representing the borders of all countries in the world. This data doesn't have to be very precise, it is to draw the countries on a computer screen. I did a bit of research and keep stumbling upon binary formats for this data to be opened in various programs and what not. Here is an example of what would be a good resource, but unusable to me.

I don't want to have to write a binary format parser if possible. So my question is two-fold:

Are there any resources for simplified country borders in plain text? (XML, JSON, CSV, etc)
If not, are there any programs where I could export the binary data found all over the web in those formats?


texts: ogr2ogr TPS spline method creates progressive error

texts:
I am trying to adjust a shapefile using ogr2ogr spline transformation with cmd commands like:

C:\OSGeo4W64\bin\ogr2ogr.exe -f "ESRI Shapefile" C:\path\output.shp -tps --optfile C:\path\gcp.txt C:\path\input.shp

I have more than 1000 control points (so they are in a separate file). And I have strange problems with precision of this method. I already saw in this question that ogr2org spline method is not really exact. But with my number of GCP and the extent of my dataset I see that the precision dramatically decreases from north to south. Like this:

In the north the method is almost exact (0.001 m error), then it smoothly loses the precision and in the south it creates an error about 60 m.

I calculated RMSE for every GCP and plotted it versus coordinates and the ID number of control point (I was creating the GCP starting mostly from north). And I have:

I tried to find and read the source code of gdal (I found gdal_tps, thinplatespline and ogr2ogr_lib modules) but I don't know that language (C++?) and don't understand how the method works. Polynomials 1, 2 and 3 order of ogr2ogr are working fine (they are not exact methods but the error is not progressing).

So, why spline precision logarithmically decreases depending on the Y coordinate? (For X coordinate I see jumps in precision every 16000 m). How is this possible? How this adjustment method works? How can I solve this problem? (I have Windows 7, 64 bit)


texts:
The process for this task so far has been to copy and paste a value from a table to the 'select by attribute' tool, then zoom to the set of selected features, make a decision about that feature/those features, update appropriate fields, and repeat. I would like to be able to run a model that will select and zoom to selected features, allow me to update some records, then automatically select and zoom to the next set of features (based on a table).

I can build a zoom to selected feature tool. I would like to shorten the copy, paste, edit time. It seems like using an iterator would be a good way to go. Is there a way to pause an iterator in ModelBuilder to allow a user to update layer attribute values - say by prompting the user with eg a 'do you want to continue?' option. If the user says yes, MB goes to the next iteration. If the user does not say yes, can MB break from the iteration loop?


texts: Overwriting “ArcGIS Online World Geocoding service” on search widget?

texts:
I am new and just exploring ArcGIS API

Wondering how to overwrite the default "ArcGIS Online World Geocoding service" each time i search it calls that service.

I wanted to call our own geocoding service. to add up. I can't also call my MapServer query using the sandbox using the search widget, but the sample 'Query' sandbox can call our MapServer.

Is there any important thing i missed like set-up on our ArcGIS server for the search functionalities to function?


texts: Performing Suitability Analysis using ArcGIS Desktop?

texts:
I'm new to GIS (using ArcGIS 10.4) and stuck on the following problem: I need a way to find ideal locations for supermarkets that will service the most people possible based on proximity (as the crow flies).

To solve this I used Location-Allocation Analysis (Network Analyst Extension) with the problem type "maximise coverage".


texts: Reading value of point (x, y) in raster and changing its value to another raster using PyQGIS?

texts:
I have two raster: raster1 with values [1,2,4,8,16,32,64,128] and raster2, same extension, with all values -9999. Starting from the coordinates of a point, I would like to read its value in raster1 and to do something like this:

If in raster1 the value of the point (x, y) = 1 and the value of the point (x-1, y) = 16, in the raster2 at point (x, y) I leave the -9999 value.

If the value of the point (x, y) = 1 in the raster1 and the value of the point (x-1, y) = 4 then in raster2 I'd write the value 1.

As output I should have the raster2 containing the values 1 only in points that interest me. Can someone help me?

This is what I have done so far:

import processing
import gdal

layer_1 = iface.addRasterLayer("path/to/raster1.tif")
if layer.isValid():
    QgsMapLayerRegistry.instance().addMapLayer(layer)

processing.runalg('gdalogr:rastercalculator', 'path/to/raster1.tif','1',None,'1',None,'1',None,'1',None,'1',None,'1','A*0-9999',None,5,None,'path/to/raster2.tif')
layer_2 = iface.addRasterLayer('path/to/raster2')
QgsMapLayerRegistry.instance().addMapLayer(layer_2)

x = 603608
y = 4359398
point = QgsPoint(x,y)
ident_1 = layer_1.dataProvider().identify(QgsPoint(x,y),QgsRaster.IdentifyFormatValue)
ident_2 = layer_2.dataProvider().identify(QgsPoint(x,y),QgsRaster.IdentifyFormatValue)
if ident.isValid():
    a = ident_1.results()
    b = ident_2.results()
    if a == 1 &&...


texts: replaceDataSource for Tables from one SDE server to another server

texts:
I found the code that works to change the layers in my mxd from one SDE workspace to the other SDE worksspace. The script runs, but does not change the tables to the new SDE workspace. I replaced lyr with Table, but can't seem to get it to work. Once I get the table data source script to update I will put it into code that loops through MXD files in a specified folder and updates any GDSERVER data connections to point to GDTEST using the assigned connection properties. I am new to this and could use some help please :)

The error msg I get is.
Traceback (most recent call last):
File "C:\Users\ferencek\Desktop\TEST1.py", line 4, in
for Table in arcpy.mapping.ListTables(mxd):
AttributeError: 'module' object has no attribute 'ListTables'

import arcpy;

mxd = arcpy.mapping.MapDocument(r"N:\\5-Processes\\Resources\\DynamicData\\SourceMXD\\Copies for GDTEST\\Test\\Data_CCTV.mxd");
for Table in arcpy.mapping.ListTables(mxd):
if Table.supports("DATASOURCE"):Table.replaceDataSource("C:\\Users\\ferencek\\AppData\\Roaming\\ESRI\\Desktop10.2\\ArcCatalog\\sde_appreader (SDE_LIVE) on GDTEST.sde", "SDE_WORKSPACE", Table.datasetName, False)

mxd.saveACopy(r"N:\\5-Processes\\Resources\\DynamicData\\SourceMXD\\Copies for GDTEST\\Test\\CopyData_CCTV.mxd");
del mxd


texts:
I drew some lines in a GIS application in this projection:
+proj=longlat+ellps=WGS84 +datum=WGS84 - is there a name for this projection?

I imported my lines into R and to prepare for calculating distance between points on my lines, I projected my lines using:
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs - again, is there a name for this projection?

To actually calculate distances between points on my lines, I used the gLength function in the R package rgeos. What sort of distance calculation does gLength use?


texts: Seeking program to highlight and calculate multiple areas on satellite image?

texts:
I would like to calculate the area of all the parking spaces in my city (North Miami Beach). I would highlight the parking in two layers, commercial and residential, and then export it to Excel and calculate the area.

What would be an easy program to do this in? Something that I can save and is easy to manipulate/correct/edit.

I am new to GIS - right now I'm looking at QGIS, but I'm not sure if this is a good program to use.

I was planning on tracing the footprints of the parking lot manually, because image recognition seems like something very complex. I could start with a small area and if the results are good and everything works out I could move to an image recognition program for a larger area.

I was planning on going over an image of google satellite maps or an equivalent.

(Now that I have been playing around with qgis I see that i don't need excel)

The idea is to calculate how much all the parking space is worth and to see how much money could be generated if we reduced driving by let say 30%. I would take the average cost of commercial and residential property and multiply it by the area.

I understand this is a rough calculation.

30% of driving trips could be replaced by cycling or reduced with the use of future car sharing technology...I would like to have a economic argument for this besides for the environmental.

(cycling would also take up space, but less then cars)


texts: How to create a list of spatial objects and merge them in R

texts:
I’m trying to load and merge a bunch of .tif spatial files from ASTER. Sorry, I can’t provide a replicable example, but I hope this ilustrates the situation:

library(raster)
pathz <- "C:/USERS/PATH/TO/WORKING/DIR/"
# The subdirectory where my .tif files are
filez <- list.files(path = paste(pathz, "/Granules", sep = ""), pattern = ".tif")
# I get the list of files
namez <- paste("r", c(1:length(filez)), sep = "")
# I create the bunch of spatial objects
for(i in 1:length(filez)){
   temp <- raster(paste(pathz, "/Granules/", filez[i], sep = ""))
   assign(namez[i], temp)
}
# This way, I've created 200 RasterLayer objects, for example
r49
class       : RasterLayer
dimensions  : 3601, 3601, 12967201  (nrow, ncol, ncell)
resolution  : 0.0002777778, 0.0002777778  (x, y)
extent      : -76.00014, -74.99986, 1.999861, 3.000139  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Paolo\Desktop\Aster DEM\Granules\ASTGTM2_N02W076_dem.tif
names       : ASTGTM2_N02W076_dem
values      : -32768, 32767  (min, max)
class(namez)
[1] "character"
namez1 <- as.list(namez)
class(namez1)
[1] "list"
namez1
[[1]]
[1] "r1"

[[2]]
[1] "r2"

[[3]]
[1] "r3"
# etc

# Obviously it did not work
m <- do.call(merge, namez1)
Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column
# I tried...
namez1 <- as.environment(namez1)
Error in list2env(list("r1", "r2", "r3", "r4", "r5", "r6", "r7", "r8",  :
  names(x) must be a character vector of the same length as x
# Also...
list2env(namez1, envir = NULL)
Error in list2env(namez1, envir = NULL) :
  names(x) must be a character vector of the same length as x
# Of course...
is.vector(namez1)
[1] TRUE

Specifically, I wish to know how to create a list of spatial objects so I can merge them.


texts: How to make new field with new characters

texts:
I want to add 0 (zero) to existing field.
I have textual field with numbers, sort ascending, 1-100:

1
2
...
13
14
...
99
100
I want to make field like this:
001
002
...
013
014
...
099
100

Is it possible to do this with some function (trim or somehow)?


texts:
I'm writing a PyQGIS script. I succeed in accessing GIS vector data (lists) from my remote server (thanks to urllib). I succeed in displaying geometries and attributes.

Now, I would like manage layer's editions and modifications. I must find a way to know the polygons that were been modified, deleted or added. Could you throw light for me?

Secondly, I would like to send these edits and these modifications toward the remote server. Could you tell me how to do that please?

In order to be more explicit I paste the useful github link:

https://github.com/picardie-nature/clicnat-qgis/blob/master/plugin.py


texts: How to ensure that Tilemill interprets my postgis query as intended?

texts:
One of my layers (a postgis query from my local osm2pgsql database) is being translated incorrectly in Tilemill.

The query is attempting to select the way, the area of the polygon, and the sport=* tag as follows:

( SELECT way, sport AS sport, way_area AS area, COALESCE(leisure, "natural", highway, amenity, tourism) AS type
  FROM planet_osm_polygon
  ORDER BY way_area DESC
) AS data

This query works as expected in pgadmin and returns the results that I expect.

In my project.mml file, this layer is represented as follows:

  {
        "Datasource": {
          "dbname": "cleveosm2pgsql",
          "extent": "-20037508.34 -20037508.34 20037508.34 20037508.34",
          "geometry_field": "way",
          "host": "localhost",
          "id": "landuse",
          "key_field": "",
          "port": "5432",
          "project": "osm-bright-imposm",
          "srs": "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over",
          "table": "( SELECT way, sport, way_area AS area, COALESCE(leisure, \"natural\", highway, amenity, tourism) AS type\n
  FROM planet_osm_polygon\n  ORDER BY way_area DESC\n) AS data",
          "type": "postgis",
          "user": "skors"
                },
        "class": "",
        "geometry": "polygon",
        "id": "landuse",
        "name": "landuse",
        "srs": "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over",
        "srs-name": "900913",
        "status": "on",
        "advanced": {}
      },

When I open my project in Tilemill, I receive the following error:

  Postgis Plugin: ERROR:  column "sport" does not exist LINE 1: SELECT
  ST_AsBinary("way") AS geom,"sport","type" FROM ( SELE...
                                            ^

  in executeQuery Full sql was: 'SELECT ST_AsBinary("way") AS
  geom,"sport","type" FROM ( SELECT way, way_area AS area,
  COALESCE(landuse, leisure, "natural", highway, amenity, tourism) AS
  type   FROM planet_osm_polygon   WHERE way_area > 10000   ORDER BY
  way_area DESC ) AS data WHERE "way" &&
  ST_SetSRID('BOX3D(-9089279.907446878
  4999593.14607681,-9030576.269723864 5058296.783799823)'::box3d, 900913)'

I have already cleared my cache (within the layer's menu in tilemill) and restarted tilemill. Any Idea why Tilemill is modifying my query to change the sport column from sport to "sport"?


texts: OpenLayers2 how to get feature info on form submit?

texts:
I have an OpenLayers pop-up with a form on selectFeature (WFS-T). Can I recover the feature.attributes in the form's onSubmit?

e.g: the form in the popup

var content = '<form onSubmit="processForm(this.form)">' +
    'Träd-Id
<input type="text" name="trädnummer" value="' +
    feature.attributes.trädnummer + '">' +
    '
Art
<input type="text" name="art" value="' +
    feature.attributes.art + '">' +
    '
Anmärk.
<textarea name="anmärkningar" rows="4" cols="40">' +
    feature.attributes.anmärkningar + '"</textarea>' +
    '
Åtgärd
<textarea name="anmärkningar" rows="4" cols="40">' +
    feature.attributes.rekommendationer + '"</textarea>' +
    '
<input type="submit" value="Spara""></form>'

and the onSubmit callback:

function processForm(oForm) {
  var feature = **recover currently selected feature here?**
  feature.attributes.trädnummer =
    oForm.elements["trädnummer"],value;
  feature.attributes.art =
    oForm.elements["art"],value;
  feature.attributes.anmärkningar =
    oForm.elements["anmärkningar"],value;
  feature.attributes.rekommendationer =
    oForm.elements["rekommendationer"],value;
  return false;
}

If so, I thought I could use WFS-T to update the attribute info? (instead of e.g. a post to PHP)
Any help appreciated


texts:
I try to use MOLUSCE plugin (versions 3.0.13, 3.0.11, 3.0.10, and 3.0.9) in QGIS 3.12.1, none of them are functioning.

Your suggestion would be of great help.


texts:
I am using the following update cursor, but when I execute it, it does not error, but also makes no changes. I am trying to Update a field where any NULL, Blank, and NEW values will be updated with 'New'.

with arcpy.da.UpdateCursor("contacts", 'Contact_Status1') as cursor:
    for row in cursor:
        if  row[0] is None or ' ' or 'NEW'  :
            row[0]= "New"
            cursor.updateRow(row)


texts:
I am trying to build plugin that allow to quick update attribute table of selected features. User select feature from vector layer then select values from plugin, click apply and table is updated. I got four combo box from which you can select values and click button "Apply changes" that should execute code to change attribute table.

My code below:

def run(self):
        """Run method that loads and starts the plugin"""

        self.dlg.connect(self.dlg.apply,SIGNAL("clicked()"),
                         self.execute_changes) #SIGNAL is not defined

        print str(result)

        if not self.pluginIsActive:
            self.pluginIsActive = True

            #print "** STARTING NM_quick_updator"

            # dockwidget may not exist if:
            #    first run of plugin
            #    removed on close (see self.onClosePlugin method)
            if self.dockwidget == None:
                # Create the dockwidget (after translation) and keep reference
                self.dockwidget = quick_updateDockWidget()

            # connect to provide cleanup on closing of dockwidget
            self.dockwidget.closingPlugin.connect(self.onClosePlugin)

            # show the dockwidget
            # TODO: fix to allow choice of dock location
            self.iface.addDockWidget(Qt.TopDockWidgetArea, self.dockwidget)
            self.dockwidget.show()

    def execute_changes(self):

            #Execute changing attributes when button is clicked
            layer = self.iface.activeLayer()
            selected_features = layer.selectedFeatures()

            fields = layer.pendingFields()

            field_names = [field.name() for field in fields]

            print field_names

            #selectedLayerIndex = self.dlg.currentIndex()
            #selectedLayer = layers[selectedLayerIndex]
            #selFeatures = selectedLayer.selectedFeatures()
            selFeatures =selected_features
            selectedLayer = layer
            ids = [f.id() for f in selFeatures]
            print ids

            selectedLayer.startEditing()

            #here i take values from combo boxes and use them to update columns
            for fid in ids:
                print fid

                selectedLayer.changeAttributeValue(fid, 6, int(self.dlg.Name_1.currentText()))
                selectedLayer.changeAttributeValue(fid, 2, str(self.dlg.Name_2.currentText()))

            print "completed changing values"
            selectedLayer.commitChanges()
            print "commiting changes"
            mc =  self.iface.mapCanvas()
            mc.refresh()
            print "refresh map"

I can access to my combo box inside changeAttribiuteValues but I would like to update table with given name not using FID number. I created list of my shapefile columns but I don't know how to use them with combo boxes.

Also I don't know how to run my code when push button is clicked. Now on when I load my plugin I got error: 'global name 'SIGNAL' is not defined'

Can someone help me?

I am using QGIS 2.16.


texts: QGIS Field Calculator Syntax Error Problem

texts:
I never before worked with the field calculator and this problem is driving me nuts ...

Interpretationserror:
syntax error, unexpected NAME
syntax error, unexpected ')', expecting $end
My expression is: ((“crimes_2020” - ("crimes_2017_2019” / 3)) / ("crimes_2017_2019” / 3)) * 100

Also, the Ok Button always seems to be greyed out.


texts: Transforming NetCDF into JSON using Python for dynamic display (with movement) using JavaScript

texts:
I have NetCDF data and am trying to transform this into JSON, using Python, to then use in a web sig (NCWMS server).
By transforming this into JSON, could I display this same data dynamically (with movement) with the use of JavaScript?


texts: Change the scene mode in Cesium every interval

texts:
I have a task which
I want to change the scene mode in Cesium every 1 mins.
For example, let say be default the scene mode is 3D. after 1 mins, it will change to 2D mode and then anohter. I used refresh feature to change the scene mode. Here is my code below.

    <!DOCTYPE html>

      <meta charset="utf-8">

      <meta http-equiv="X-UA-Compatible" content="IE=edge">

      <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no">
      <title>Hello World!</title>

It does not work except my console.log.

Can you help me?


texts: TopologyException: side location conflict?

texts:
I'm receiving the following error on a query:

  ERROR: GEOSIntersects: TopologyException: side location conflict at -121.12159399072365 36.206327556592377

I've tried running ST_isValidReason() as suggested in the accepted answer on a related question here but all of the geometry fields are returned as valid. I don't have any null geometry either.

How else can I troubleshoot what may be causing this error?


texts: Purpose of geography_columns and geometry_columns views in PostGIS?

texts:
Can anyone tell me what is the purpose of geography_columns view and geometry_columns view in a PostGIS database?

My geography_columns view is empty so I can't see what the purpose of it is.


texts: How to implement the 'contain' spatial join in Google Earth Engine?

texts:
How to select polygons which are contained in another polygon in Google Earth Engine? For example, in this case below, how to select only the red square completely within the green square, excluding the red square intersected?

As far as I know, the functions, .filterBounds() or the combination of ee.Join.saveAll() and ee.Filter.withinDistance() , are utilized to conduct the spatial query of 'Intersection', not 'Contained'.


texts: Deleting blank fields in multiple shapefiles in folder using ArcPy?

texts:
I have some shapefiles and each shapefile is having some blank or zero value fields. My task is to delete blank/zero value  fields for all shapefiles in the folder. I have searched for solution in this site and found one Python program but unfortunately I am not successful with available Python program. Below is the my code and I don't know where I made a mistake.

I have attached image also for reference and the fields i wanted to delete are shown with line

After running code, I did not encounter any error message but fields are not removed

result:

Empty field deleted

import os
import arcpy
from arcpy import env
path = r"C:\SOFTWARE\Shapefiles\test"
env.workspace = path
fcs = arcpy.ListFeatureClasses()

# Iterate over the feature classes
for fc in fcs:
    # Test whether the field is required--rather than testing the field type
    fieldList = arcpy.arcpy.ListFields(fc)

    # Iterate over the fields
    for field in fieldList:
        # Set a test variable to False
        with arcpy.da.SearchCursor(fc, field.name) as cursor:
            for row in cursor:`enter code here`
                if field == " " or field == 0:
                    arcpy.DeleteField_management(fc, field.name)
print ("Empty field deleted")


texts: Exposing GRASS GIS add-on in QGIS Processing framework?

texts: Having QGIS and GRASS GIS 7 installed, how can I expose a GRASS GIS add-on, which has been installed using g.extension, in QGIS' Processing framework?


texts: ST_Within SQL returning too many results

texts:
This ST_Within query is working fine

SELECT school.school_name, Courses.Course_Title, Students.No_of_Students, district.name
FROM (school INNER JOIN Students ON school.school_id = Students.School_id) INNER JOIN Courses ON Students.Course_id = Courses.Course_id, district
WHERE district.name = 'Rochdale'
AND ST_WITHIN(school.geom,district.the_geom)

But when I try and return another district along with 'Rochdale' it just returns loads of results

this is what I am using

SELECT school.school_name, Courses.Course_Title, Students.No_of_Students, district.name
FROM (school INNER JOIN Students ON school.school_id = Students.School_id) INNER JOIN Courses ON Students.Course_id = Courses.Course_id, district
WHERE district.name = 'Rochdale' or district.name = 'Calderdale'
AND ST_WITHIN(school.geom,district.the_geom)

Can anyone help?


texts: How to add complex Emergency Symbology (TrueType fonts) to QGIS?

texts:
In other packages I would simply ensure the TrueType font (+.ttf) was loaded on the operating system and it would be available for use.  QGIS seems to rely on making your own. The symbols I wish to import are complex Emergency Symbology that couldn't be recreated within QGIS.

How do I add complex Emergency Symbology (TrueType fonts) to QGIS?


texts: Get the geotransform in Mercator given a pixel size away from the origin

texts:
I have a radar image with the following information:

The projection is Mercator (Mercator_1SP in GDAL)
The projection origin is 0º lat 0º lon
At a given geographic coordinates (e.g. lon=10º lat=40º), the scale is 1kmx1km each pixel. I know at wich pixel in the image this coordinates are (they are usualy at the center of the image)
The datum can be WGS84 or others, so it's not always the spherical earth

How can I determine the geotransform? I don't know how to transform this pixel scale distance into the projection's scale at the given point.

--EDIT--

So I need the upper left pixel coordinates and the pixel width and height in the Mercator coordinates from the data explained above.

Thank you.


texts: I am creating an interactive map using ArcGIS online of which I will share the link for people to view and interact with the map. I have 7 layers and will want the viewers to be able to turn on as well as turn off layers as they please. However, when I click on the generated link to test it, the layers appear fixed and I can`t turn it on and off. How do I go about this?


texts: XML Filter Expression for WFS download using Interoperability Connection Extension in ArcGIS Desktop

texts:
I need to download features from this WFS but there are too many features for the server to handle my request. So I thought I should filter my request using the XML Filter Expression Editor. However I cannot find any information on how to write such an expression into the editor.

I only need features with the attribute "FBL" having the value '8' (it is a string and not a numerical field).


texts: What technology does google analytics use to show the map?

texts:
Does anyone know what technology is used by google analytics in the Geo location tab?

It looks like a combination of image and vector. Is it WMS & WFS or something else?


texts: Placement of Save Button in ArcGIS Javascript AttributeInspector dijit

texts:
I am having a problem with an AttributeInspector that I noticed also exists in an online sample here: http://developers.arcgis.com/javascript/samples/ed_attribute_inspector/

If you click on a polygon to display the AttibuteInspector, you will notice that hovering your mouse over the upper right of the delete button selects the save button. In the Chrome console, I was able to find that there is a rogue  tag that might be causing this issue. See attached screenshot.

The code that created this is as follows and is seen in the sample code here.

var saveButton = new Button({ label: "Save", "class": "saveButton"});
domConstruct.place(saveButton.domNode, attInspector.deleteBtn.domNode, "after");

How can this issue be avoided so that no parts of the two buttons or their containers overlap?


texts: Interpret terrain type from height data and/or geolocation

texts:
I have downloaded a large subset of the STRMGL1 Dataset, and have further split each entry down into 5km by 5km tiles. I now need to start classifying each tile by terrain types. Are there any automated methods for doing so from the height data and/or the Geolocation of each tile? For example, is there a database of terrain types at different locations?

Any terrain types will work, but I would prefer to have at least some of the following:
OPEN - TUNDRA - OASIS - STEPPE - HILL - DESERT - MOUNTAIN - FOREST - MARSH - BOG - SWAMP

Multiple tags per tile are allowed.


texts: Sharing auxiliary data between projects using QGIS?

texts:
I am using QGIS 3.2.0 and I want to share auxiliary data between projects but I could't find how to do that.

I tried to "play" with auxiliary layer export but there is no auxiliary layer import.

Also I tried and just duplicated and renamed .qgd file of the first project but no solution.

See screenshots:


texts:
I have problems to visualize ecw files for qgis 2.3.0 using ubuntu 12.04 LTS. From terminal, following the instructions for qgis 1.8 found at ECW for QGIS 1.8 in Ubuntu 12.04, because launching the command sudo gdal-ecw-build /usr/local/.
I added the repository ppa:ubuntugis/ubuntugis-unstable and then did these steps:

wget http://de-mirror.org/gentoo/distfiles/libecwj2-3.3-2006-09-06.zip
unzip libecwj2-3.3-2006-09-06.zip
cd libecwj2-3.3
./configure
make
sudo make install
sudo gdal-ecw-build /usr/local/

Launching the last command I receive some warnings like this:

ecwdataset.cpp:1501:15: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]

and then an error message:

/usr/local/include/gdal_priv.h:225:19: error: ‘char**
GDALOpenInfo::papszSiblingFiles’ is private
ecwdataset.cpp:2526:48: error: within this context

Finally (in Italian because I have an Italian version of Ubuntu 12.04 LTS):

make: *** [ecwdataset.o] Errore 1

Can anybody help me?


texts:
I read a lot of speeding up the browser widget, but nothing seemed to help me.
as long i have it closed qgis is fast but if I open it all actions slows down.
My project is on a local ssd - fast dir scan is enabled on all devices.
what can I do?
regards.


texts: Set Up WFS Service on MapServer and display in OpenLayers

texts:
I have used following map file to create the WFS service on MapServer but struggling to display the service output using OpenLayers.

Please suggest some links to use WFS service in OpenLayers.

MAP
  NAME "WFS_server"
  STATUS ON
  SIZE 600 300
  UNITS DD
  IMAGECOLOR 255 255 255
  #
  # Start of web interface definition
  #
  WEB
    IMAGEPATH "/ms4w/tmp/ms_tmp/"
    IMAGEURL "/ms_tmp/"
    METADATA
      "wfs_title"          "WFS Demo Server for MapServer" ## REQUIRED
      "wfs_onlineresource" "http://demo.mapserver.org/cgi-bin/wfs?" ## Recommended
      "wfs_srs"            "EPSG:4326 EPSG:4269 EPSG:3978 EPSG:3857" ## Recommended
      "wfs_abstract"       "This text describes my WFS service." ## Recommended
      "wfs_enable_request" "*"  # necessary
    END
  END
  PROJECTION
    "init=epsg:4326"
  END

  #
  # Start of layer definitions
  #

  ##################
  # World Continents
  ##################
  LAYER
    NAME "Mumbai_City"
    METADATA
      "wfs_title"         "World continents" ##REQUIRED
      "wfs_srs"           "EPSG:4326" ## REQUIRED
      "gml_include_items" "all" ## Optional (serves all attributes for layer)
      "gml_featureid"     "ID" ## REQUIRED
      "wfs_enable_request" "*"
    END
    TYPE POLYGON
    STATUS ON
    DATA '/ms4w/data/Mumbai_City_Boundry.shp'
    PROJECTION
      "init=epsg:4326"
    END
    CLASS
      NAME 'Mumbai City'
      STYLE
        COLOR 255 128 128
        OUTLINECOLOR 96 96 96
      END
    END
  END #layer
END #mapfile


texts:
I have many polygons with overlaps.

I want to remove the part of overlapping that is between polygons.

I tried Symmetrical difference, it gives me only the area overlapping and not two polygons.

Red polygons are overlapping


texts: Working with 2D building footprints and ArcGIS Pro / CityEngine

texts:
I've been trying for a few days now to create 3D buildings from my 2D building footprints. All the tutorials I've been through let me down right where it counts so I've decided to ask the pro's here.

My 2D footprints do have height attributes and I know there is a way to extrude them. However I cannot make this happen with a single building nor the entire group. Can someone explain in fine detail how I should get this done?

The image below is after selecting one building. You can see that height exists, however extruding rules are not explained very well here.


texts: Install php mapscript on ubuntu 12.04 PHP 5.4

texts:
I follow this website : http://blog.geoghana.com/2013/06/28/setting-up-lampp-linux-apache-mapserver-php-postgresqlpostgis/ to install mapserver and the php extension on my ubuntu 12.04 64bit and PHP 5.4 system (i use this ppa https://launchpad.net/~ondrej/+archive/ubuntu/php5-oldstable). I almost succeed until i have a problem on installing php5-mapscript :

sudo apt-get install php5-mapscript

Here is the message on my linux terminal:

  Some packages could not be installed. This may mean that you have
  requested an impossible situation or if you are using the unstable
  distribution that some required packages have not yet been created
  or been moved out of Incoming.
  The following information may help to resolve the situation:

  The following packages have unmet dependencies:
   php5-mapscript : Depends: phpapi-20090626


texts: Image disappearing when georeferencing google earth image?

texts:
I have downloaded an image from google earth with four control points ready to geo reference. I enter the first two points and everything is fine. Then I enter a third point, the image disappears and the 'zoom to layer' just shows the control points with a blank map behind them.

Has anyone else experienced this and would know how to rectify it? I don't know if its a problem with the control points or the image itself.


texts: Clipping SQL Server LINESTRING geometry based on polygon without overlap in result?

texts:
Using DbGeometry and Entity Framework to work with Sql Geometry types.
User can draw a line that can run outside of a given polygon. my validation traps that. I correct this by clipping the line to be the intersection of both geometries.

And that is working. however, the result of that intersection still trips up my validation (i assume by one coincident point).

I either need to shorten that line by a tiny amount or redesign my validation rule(which uses Shape.Within() to determine if the line is outside the polygon).
I have solved this issue with polygons in the past by buffering them -0.00001 but that doesn't appear to be an option for lines


texts:
I have geoserver layer with sld style applied. The sld have like 50 rules using type attribute of the feature. Now i need to generate legend but only for specific area (i have bbox coordinates). I am doing this on server side and i can get all the features and their attributes of the choosen area (but i don't have access to style graphics ).

The only way i can think of implementing this is to generate sld style (rules) dynamically and pass it to geoserver as SLD_BODY. It seems a lot of work so I wonder maybe there is a better way of doing this?

Or maybe it is possible to pass some kind of parameter array to getLegendGraphic request like array of features ids or their type?


texts: Calculating true distance of cost distance raster using ArcMap

texts: I have calculated cost distance using a raster unrelated to path distance (a heat index) but would like to also measure the true path distance of those 'routes.' Is there a way to calculate such? If I simply calculate path distance, the tool obviously takes different paths than that used for cost distance. I'm using ArcMap.


texts: Cannot enable offline mode for ArcGIS Online web map in order to use in Esri Collector app

texts:
I've been trying to enable offline mode in a WebMap in order to use one Collector. I reviewed the requirements here however I enabled sync and caching on my feature layer and map layer I believe. Here is the option in ArcServer Manager that shows it as enabled:

What else can I check to ensure that I meet all the requirements in order to enable offline mode in WebMap? And how can I accomplish this. I see all the tutorials online refer to ArcMap, however, I use ArcGIS Pro and I'm having a hard time finding a tutorial that matches the products I'm using.

EDIT: I used the World Imagery Basemap for export from ArcGIS online from the tut included in the link above. This is th link to the basemap for Export: arcgis.com/home/item.html?id=226d23f076da478bba4589e7eae95952

Also I should mention that I want my data to upload to a geo-database when bringing in the data into the office for syncing.


texts: How to infer a relationship between two spatial points using R

texts:
I am quite new to the area of spatial statistics, but I'm very interested. For learning and demo purposes, I've created three datsets.
Dataset - Persons: This describes individuals at a certain location with a few variables. Please note, that the persons are located in the provided cities. A short explanation:

POINT_X: X-coordinate of city.
POINT_Y: Y-coordinate of city.
city: The name of the city, where they live.
ill: "1"  states that they are ill. For learning purposes, all persons are ill.
job: If they have a job or not. "1" means: they have one, "0" means they haven't got one.
disnw: The distance to the nearest waterpoint.
wID: not relevant.

Dataset - City: This describes a number of cities including some variables. A short explanation of these:

city: The name of the city.
population: The population of the city.
POINT_X: X-coordinate of city.
POINT_Y: Y-coordinate of city.
ill: Number of ill persons in the city.
notill: Number of healthy persons in the city.
disnw: The distance (in km) to the nearest waterfeature.
wID: not relevant
rate_ill: The rate of ill persons in the city.
rate_notill: The rate of healthy persons in the city.

Dataset - Waterfeatures: . Please note that the viallages are on the same location as persons. This is a collection of spatial points, which describes waterfeatures.

POINT_X: X-coordinate of a waterfeature.
POINT_Y: Y-coordinate of a waterfeature.

geographic overview about the setting (red are persons, blue are waterfeatures, yellow are cities)

Now I want to check the hypothesis that cities, which are nearer to waterfeatures (so where the variable disnw is lower), have a higher number of ill persons. So is there a correlation between the number of ill persons(no_ill)/rate of ill (rate_ill) and the proximity to water features (disnw). I know, that the datasets are possibly not representative or suitable for my hyptothesis, but for now this fact shouldn't matter.
I've already looked at some functions and packages, but I'm very unsure about a suitable method. Methods, which might be useful (at least from my point of view): semivariogram, variogram, Ripley's K function, G-Function, correlation coefficient.

To give you a better overview, I've created example datasets. You can find these here:

persons = read.csv("http://pastebin.com/raw.php?i=3aMGi9Ax", header = TRUE, stringsAsFactors=FALSE)
city = read.csv("http://pastebin.com/raw.php?i=Lk3KXLQT", header = TRUE, stringsAsFactors=FALSE)
water = read.csv("http://pastebin.com/raw.php?i=hQRvMZwE", header = TRUE, stringsAsFactors=FALSE)

It would be awesome to get some input from your side. Maybe you have a tip, how to perform this kind of analysis.
Thanks in advance!


texts: Comparing NDWI regions with a given shapefile Google Earth Engine Python?

texts:
With the code below, I computed the NDVI, NDWI and NDWIMOD indices of a R.O.I. in order to detect water areas in that R.O.I.

Now, my goal concerns to compare the water areas with a given shapefile in order to find which areas are correctly detected by computing the indices. I want to do it for each index separately.

Do you have any advice about how I can carry out it? Moreover, How can I save the new areas detected?

The code is the following:

import ee
import ee.mapclient

ee.Initialize()
ee.mapclient.centerMap(15.88, 41.78, 10)

def maskS2clouds(collection):
    qa = collection.select('QA60')
    cloudBitMask = 1 << 10
    cirrusBitMask = 1 << 11
    mask = qa.bitwiseAnd(cloudBitMask).eq(0)and(qa.bitwiseAnd(cirrusBitMask).eq(0))
    return collection.updateMask(mask).divide(10000)

polygon = ee.Geometry.Polygon([[
[16.440171600655958,42.033470726105556],
[15.511826874093458,42.12316884368949],
[14.621934295968458,41.689821328579825],
[15.441789032296583,41.00320702735848],
[16.246537567452833,41.36906568596335],
[16.440171600655958,42.033470726105556]]])

collection = ee.ImageCollection('COPERNICUS/S2').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)).filterDate('2018-01-1','2018-06-30').map(maskS2clouds).filterBounds(polygon)

medianpixel = collection.median()
medianpixelclipped = medianpixel.clip(polygon)

ndvi = collection.mean().normalizedDifference(['B8','B4'])
ndviViz = {'min':0.6, 'max': 1, 'palette': ['0000FF', 'FF0000']}
ndvimasked = ndvi.updateMask(ndvi.lt(0.0))
ee.mapclient.addToMap(ndvimasked, ndviViz, 'Sentinel_2 mosaic')

ndwi2 =  collection.mean().normalizedDifference(['B3','B8'])
ndwi2Viz = {'min': -0.4, 'max': 0.0, 'palette': ['FF0000', '00FF00']}
ndwi2masked = ndwi2.updateMask(ndwi2.gt(-0.125))
ee.mapclient.addToMap(ndwi2masked, ndwi2Viz, 'Sentinel_2 mosaic')

ndwi2mod = collection.mean().normalizedDifference(['B8','B3'])
ndwi2modViz = {'min': 0.0, 'max': 0.4, 'palette': ['FFFFFF', 'FF0000']}
ndwi2modmasked = ndwi2mod.updateMask(ndwi2mod.lt(0.0))
ee.mapclient.addToMap(ndwi2modmasked, ndwi2modViz, 'Sentinel_2 mosaic')

capitanata = ee.FeatureCollection('users/vpaparella96/vasche_puglia')
ee.mapclient.addToMap(capitanata)

I shared to everyone the shapefile.


texts: ArcGIS javascript API Printing: fuzzy result

texts:
I created a real simple web app, copied from samples actually, to try figuring out something about the Print Dijit.
Here's the thing. When you click print here there are two major concers with the result:

The symbology, especially the rectangle in the bottom left corner, are mashed up.
The background layer is not that high quality as it is in the browser

Any thoughts how this is possible and how to resolve both issues?


texts: Error when creating Polyline from parsed table values: "CreateObject cannot create geometry from inputs"

texts:
I'm writing a script that parses GPS events from an SDE table to try and create lines in a feature class representing when street sweeper brushes were turned on.

I'm appending the lat and long into respective lists then trying to use an insert cursor to create the lines.

cursor = arcpy.da.SearchCursor(tbl, ["TIME", "EVENT", "LATITUDE", "LONGITUDE"])
icursor = arcpy.da.InsertCursor(SweptRoadsFC, ["SHAPE@"])

XPointList = []
YPointList = []
Sweeping = 0
array = arcpy.Array()

for row in cursor:
    xpnt = float(format(row[3]))
    ypnt = float(format(row[2]))
    if Sweeping == 0:
        if format(row[1]) == "Brushes On":
            Sweeping = 1
            XPointList.append(xpnt)
            YPointList.append(ypnt)
    else:
        if format(row[1]) != "Brushes Off":
            print "Brushes are on - " + format(row[1])
            XPointList.append(xpnt)
            YPointList.append(ypnt)
        else:
            Sweeping = 0
            print "Brushes are now off - " + format(row[1])
            point = arcpy.Point()
            count1 = 0
            for x in XPointList:
                point.X = x
                point.Y = YPointList[count1]
                count1 += 1
                array.append(point)
            geometry = arcpy.Polyline(array)
            icursor.insertRow([geometry])
            array = []
            XPointList = []
            YPointList = []

Currently, I'm getting "RuntimeError: Object: CreateObject cannot create geometry from inputs" on the insertRow line.

I believe it's some sort of data type or formatting error, as when I hard-code the first three points into lists it runs fine, like this:

XPointList = [-97.930334, -97.9302865, -97.9302068]
YPointList = [29.8852428, 29.8852505, 29.8851645]

icursor = arcpy.da.InsertCursor(SweptRoadsFC, ["SHAPE@"])
array = arcpy.Array()
point = arcpy.Point()
count1 = 0
for x in XPointList:
    point.X = x
    point.Y = YPointList[count1]
    count1 += 1
    array.append(point)
geometry = arcpy.Polyline(array)
icursor.insertRow([geometry])

Any ideas on what would be throwing the error?

I've tested the conditionals with print statements and the array contains the correct coordinates when the first polyline is trying to be created.


texts: How to use ST_Difference in case of an featureless input table?

texts:
I use the query posted below to create the difference of two polygon layers.

CREATE VIEW public.v_diff AS
SELECT
  row_number() over() AS gid,
  subquery.geom
FROM
(SELECT
  (ST_Dump(COALESCE(ST_Difference(layer_a.geom, ST_Union(layer_b.geom)), layer_a.geom))).geom::geometry(Polygon, /*SRID*/) AS geom
  FROM public.layer_a, public.layer_b
  GROUP BY layer_a.gid
) AS subquery;

Some of my datasets contain no features for 'layer_b'. In this case the query returns no results. Is there any way to return the features of 'layer_a' instead?


texts: OpenLayers: How do I remove a layer based on a value?

texts:
In React, I'm using a pagination that returns a certain amount of results on each page. Using the address of these results, I convert them to longitude and latitude using a geo-coder package and set my markers.

setMarkers(map){
    const geoCode = new GeoCode();
    const dataLength = this.props.getData.length;
    for(var i = 0; i < dataLength; i++) {
        geoCode.geolookup(this.props.getData[i].address)
        .then(result => {
            // ADD A LAYER EACH TIME AN ADDRESS IS CONVERTED TO
            // LONGITUDE AND LATITUDE
            map.addLayer(new VectorLayer({
                source: new VectorSource({
                    features: [new Feature({
                        geometry: new Point(fromLonLat([result[0].lng, result[0].lat])),
                        name: 'Marker',
                        id: 'this.props.getData[i].id_number'
                    })],
                }),
                name: 'Marker', // THE NAME OF THE VECTOR LAYER
                zIndex: 5
            }))
        })
    }
    this.mapClickHandler(map)
}

However, each time the user goes to a new page, all the marker layers should be removed so that I can add new marker layers based on the new results.
I know that I have to remove layers like this:
map.removeLayers()
 But how do I use this to remove layers based on the layer name value 'Marker'


texts: Missing features from featureCollection in Google Earth Engine but present in Fusion Table

texts:
I am using google earth engine code editor to extract some band values at specific point locations. To do so I have loaded my .csv file with point IDs and coordinates into Fusion Tables and the Map looks ok, it displays all the 48455 points.

I have then imported the fusion table as a Feature Collection in Earth Engine code editor, which seems to be ok as it says that the collection contains 48455 rows. However, the points located in South East Asia and Oceania are not displayed in the map.

It seems not to be just mapping issue because when I run ReduceRegions to get the band values at the points I get NAs for those located south east asia and Australia. I have double checked the correctness of the coordinates by plotting the points in R and it worked fine. Thus, I guess it's not a projection issue either (they are all WGS84).

Link to fusion table:
https://fusiontables.google.com/DataSource?docid=1n7yJSCuigXu6J-E5_NsUwHemVPyQ9R-NeNEFe-OA#map:id=3

Link to earth engine code: https://code.earthengine.google.com/95c044ed96d7ef0a1267ed5529736ef2

Someone can help me?


texts: Conversion of GRIB and NetCDF to my database

texts:
.

I have downloaded "High Resolution Initial Conditions" climate forecast data for one day, it was in extension .tar.gz so I extracted it in my local directory and I get the files like in the attached image.

I think, that the files without extension are GRIB data (because first word in them is "GRIB"). So i want to get data from the big files (GRIB and NetCDF formats containing climate data like temerature & pressure in grid) to my database, but they are binary.
Can you recommend me some easy way for getting data from these files? I can't get any information about handling their datasets on their website.

Converting these files to .csv would be nice, but I can't find any working program to convert the GRIB files.


texts: Executing a QGIS model per feature of a polygon

texts:
I have just got into QGIS models and want to run the attached model per feature in a polygon shapefile. I have been trying to do this in two ways and neither of them work. I tried find unique---select the same value---use it to clip but I get an error.

I just want to it to run per feature in the guidelines file and then merge all the final outputs. I also want each file to be the row value. In ArcGIS models it would just be "%FieldName%_SEQ.shp" for example.

Happy to share the test data if anyone can help.


texts:
I have multiple mxd's with many company logos as jpg's.  We've moved the location of the logos from the F:\ drive to the C:\ drive.  I have a script that is supposed to change the path, but not sure where I enter the find/replace.  I did not write this script.

##Script used to change all JPEG images with old file paths to the new##

#import modules
import arcpy
import sys

#set environment
mxd = arcpy.mapping.MapDocument("CURRENT")
clean = 0

#loop through list and change all picture element file paths
for pic in arcpy.mapping.ListLayoutElements(mxd,"PICTURE_ELEMENT"):
    path = pic.sourceImage
    pathNEW = path
    pathNEW = pathNEW.replace(r"F:\Logos\", "C:\GIS\Logos\")
    try:
        pic.sourceImage = pathNEW
    except:
        clean = 1
        pic.sourceImage = path

#refresh the active view
arcpy.RefreshActiveView()

print "Auto changes complete!"
if clean == 1:
    print "One or more images need manual correction."


texts:
I am trying to generate 3d Terrain.
I am serving the terrain files at localhost:8081

Here is the Cesium Code

var viewer = new Cesium.Viewer('cesiumContainer', {
    baseLayerPicker : false,
    terrainProvider : new Cesium.CesiumTerrainProvider({
        url : 'localhost:8081'
    })
});

I get the following error:

  An error occurred in "CesiumTerrainProvider": An error occurred while accessing localhost:8081/layer.json.

This is my layer.json file and i have given appropriate access permissions:

{
  "tilejson": "2.1.0",
  "format": "heightmap-1.0",
  "version": "1.0.0",
  "scheme": "tms",
  "tiles": ["{z}/{x}/{y}.terrain?v={version}"]
}


texts: Find closest NON interpolated point on a LINESTRING

texts:
I have a linestring and a point, I want to find the closest point on the linestring to that point.

Below is a screenshot of a test I created to confirm that the various POSTGIS functions for finding the closest point between two objects will interpolate a point that isn't in the original linestring.

I need to return the nearest point that is an actually point that exists in the linestring source. In the test it should be the first square-red marker i want returned. But POSTGIS returns a point in-between points 1 & 2 that it generates.

[test screenshot]
http://i.imgur.com/XWfJPZV.jpg


texts: MODIS data with .hdf extension not compiling in Ubuntu

texts:
I am trying to use GDAL to read .hdf file, but it does not seem to recognize the extension. It works fine in Windows system.

The GDAL open command gives the following error message:

  ERROR 4: `MYD06_L2.A2014002.0645.006.2014269075802.hdf' not recognised
  as a supported file format.

You can access the .hdf file from this link [~80 MB].


texts: Spyder 3.2.6 does not work with GRASS GIS 7.2.2

texts:
I've installed spyder 3.2.6 with pip from the OSGeo console. There is a spyder.exe file in ~\OSGeo4W64\apps\Python27\Scripts and the Spyder package directories show up in site-packages. If I run the executable, the IDE does not appear and there are no error outputs in the console. Visibly, absolutely nothing happens.

I've looked at

Using Spyder IDE to develop python for Grass GIS, get 'ImportError: No module named grass.script'
How to control QGIS from within Anaconda's Spyder (Continuum Analytics Python distribution) in Windows

Following those exchanges, I do notice that I don't have a spyder.bat file, which I assumes just calls spyder.exe and some other things.

Has anyone else tried Spyder with the latest GRASS and had any success?


texts:
Came across this "Create layer with selected" button in the result for a Spatial Query in QGIS. And simply, this creates another layer from the same source with a query on the IDs. Does anyone know if this is available elsewhere in the GUI too? I looked around and couldn't find it, but it'd sure be handy.

There's been some confusion as to what I'm asking. I will try to clarify.
If there is a subset of features I want to work with, I begin with selecting them; this could be a spatial query, manual selection, or query by attributes. Once these are selected, if I don't want to create a new shapefile (or whatever), but instead work with and maybe modify the original data, I duplicate the layer and do a query on it. If it's a subset based on a common attribute, this is fairly straightforward, but if the attributes aren't common, I either have to type in every ID into the query or export the selection into a table, copy the ID column, paste it into notepad, find and replace characters to get the search string I'm looking for, copy that and paste it into the definition query with the additional search verbiage. This button does all that for me in a second. I'm not looking to create a new shapefile. I just want a defined subset to work with in the original data. That what this button does. It takes the IDs of the selected features and creates a new layer from the original source and defines a query on the IDs for the layer.


texts: Error messages when converting feature to raster

texts:
Have tried converting feature to raster multiple times but consistently getting 2 error messages:

Error 000871: Output raster: Unable to delete the output D:\xxxxxxx

Error 000581: Invalid parameters

It begins the conversion, and I get 3 .adf files, instead of 5 .adf files, a log and metadata, then outside the folder .aux and .ovr files.

I have another project in which I successfully converted to raster, but this isn't open so I don't think the error is caused by a LOCK file.

Not sure what the invalid parameters error is caused by - have done exactly the same thing before without it giving me errors. Have attached a screenshot in case it helps.


texts: Import OSM data in postgis with python script

texts: I want to import OSM data (way and nodes) in postgis database with a script pyhton, I found in the internet that I should create a schema but I don't know how and I really don't know what to put on the script and which library should I use.


texts: Installing headless ArcGIS Runtime SDK Qt Linux64 1026 on Linux Mint 19.1?

texts:
I want to install ArcGIS runtime SDK on Linux with commandline.
I can't use wizard installer.

Do you have any solution?


texts: Polygon object giving the wrong coordinates for a bounding box in ArcPy

texts:
I am working on a script where I identify polygons that contain other polygons, and I've run into a very perplexing issue. In the image below, green polygons have been identified as "container" polygons, and yellow polygons are supposedly polygons that are either not containers or are not contained.

The red highlighted polygon has holes where the smaller "polygon-itos" (see turquoise polygon) are located, but it should still be green. When checking how the script might have gone awry, I found that the polygon object was not reporting the bounding box of the container polygon correctly. I obtained maximum coordinates for each polygon using a search cursor:

polygon_dict = {}
with arcpy.da.SearchCursor(polygon_fc, ['SHAPE@', ID]) as cur:
    for row in cur:
        xmax = row[0].extent.XMax; ymax = row[0].extent.YMax
        xmin = row[0].extent.XMin; ymin = row[0].extent.YMin
        polygon_dict[row[1]] = [xmax, ymax, xmin, ymin]

polygon_dict[red_polygon] = [3072744.5390000045,
                             1254989.898000002,
                             3072732.677000001,
                             1254981.4430000037]
polygon_dict[turquoise_polygon] = [3072811.898000002,
                                   1255530.0219999999,
                                   3072794.446999997,
                                   1255488.5189999938]

turq = polygon_dict[turquoise_polygon]
red = polygon_dict[red_polygon]

When I compared the maximum coordinates of the red polygon and the turquoise polygon using turq[0] < red[0] and turq[1] < red[1] and turq[2] > red[2] and turq[3] > red[3], I expected to see that turquoise would be completely within the red, but that was not the case. Just looking at the picture indicated that the turquoise polygon should at the very least touch the red one, but even testing that theory with the polygon geometry method red_poly.touches(turquoise_poly) returned false.

Why is this happening when the turquoise polygon's bounding box is clearly completely inside the red polygon's bounding box?


texts:
I recently installed Xray for Catalog and my Catalog has been insanely slow ever since. I've uninstalled the Xray extension but my Catalog is still slow.

Any suggestions to speed Catalog back up to normal?


texts: How to create an Azimuthal Equidistant map for long-distance radio?

texts:
I'm trying to produce world or continent scale maps that use the Azimuthal Equidistant projection. I'm a radio amateur, and to an antenna, the whole world looks like spread out on a flat plane, with the antenna at the centre and the whole world at an (r, Θ) from that point.

There is software to do this - AZ_PROJ - but it uses its own database, and is entirely written in PostScript, so it's not exactly hackable with the usual geo tools.

OGR has the 'aeqd' projection, which can be invoked like this:

ogr2ogr -t_srs "+proj=aeqd  +lat_0=43.7308 +lon_0=-79.2647" out.shp in.shp

for a centre at 43.7308 N, 79.2647 W. Unfortunately, this wraps everything beyond the current hemisphere "behind" the map. It also projects the coordinates into cartesian metres, where I'd prefer polar.

Any suggestions for alternatives, please?


texts:
I have many GPX files and i want to show them as layer on mapbox tiles.
I managed to create "simple polyline" from gpx and then build GeoJson which include many features.
That GeoJson is pretty big, many trips are contain with each other.
Is there an efficient way to create one shape out of many GPX files ?


texts:
I have 6 point .shp files wich contain elevation data, the first 3 files are rather small (about 200MB for the shp and 1GB for the dbf) and the last 3 much larger (500MB .shp and 3.5GB for dbf). I would like to create filled polygons between the lines. The plugin "generate contours" is perfect for my needs.
Unfortunately, the tool only works for small files. When running the script with larger files I get the following error message.

QGIS-Version: 3.10.12-A Coruña QGIS-Codeversion: 75c848ffb1 Qt-Version: 5.11.2 GDAL-Version: 3.1.4 GEOS-Version: 3.8.1-CAPI-1.13.3
PROJ-Version: Rel. 6.3.2, May 1st, 2020 Verarbeite Algorithmus…
Algorithmus Generate Contours startet… Eingabeparameter: {
'ContourInterval' : 2, 'ContourLevels' : '', 'ContourMethod' : 3,
'ContourType' : 1, 'DuplicatePointTolerance' : 0, 'ExtendOption' : 0,
'InputField' : '"zfield"', 'InputLayer' :
'C:/Users/denni/Desktop/462000_5520000_dgm1.shp',
'LabelDecimalPlaces' : -1, 'LabelTrimZeros' : False, 'LabelUnits' :
'', 'MaxContourValue' : None, 'MinContourValue' : None, 'NContour' :
2000, 'OutputLayer' : 'TEMPORARY_OUTPUT' }
Triangulating 17936620 points Traceback (most recent call last): File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGenerator.py",
line 585, in filledContourFeatures trig,z=self.trigContourData() File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGenerator.py",
line 472, in trigContourData trig=self.buildTriangulation(x,y) File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGenerator.py",
line 462, in buildTriangulation trig=Triangulation(x,y) File
"C:\PROGRA~1\QGIS3~1.10\apps\Python37\lib\site-packages\matplotlib\tri\triangulation.py", line 54, in init self.triangles, self._neighbors =
_qhull.delaunay(x, y) RuntimeError: Error in qhull Delaunay triangulation calculation: insufficient memory (exitcode=4); use
python verbose option (-v) to see original qhull error.
During handling of the above exception, another exception occurred:
Traceback (most recent call last): File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGeneratorAlgorithm.py",
line 326, in processAlgorithm for feature in
generator.contourFeatures(): File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGenerator.py",
line 588, in filledContourFeatures raise
ContourGenerationError.fromException(sys.exc_info())
contour.ContourGenerator.ContourGenerationError: ['RuntimeError: Error
in qhull Delaunay triangulation calculation: insufficient memory
(exitcode=4); use python verbose option (-v) to see original qhull
error.\n']
During handling of the above exception, another exception occurred:
Traceback (most recent call last): File
"C:/Users/denni/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\contour\ContourGeneratorAlgorithm.py",
line 330, in processAlgorithm feedback.reportError(ex.message())
TypeError: QgsProcessingFeedback.reportError(): argument 1 has
unexpected type 'list'
Ausführung nach 23894.76 Sekunden gescheitert
Lade Ergebnis Layer Algorithmus 'Generate Contours' beendet

It says "insufficient memory" but I still have 71GB on my hard drive. It should be enough shouldn't it. I do have "only" 8GB RAM, perhaps that's not enough. Does anyone have any ideas on what I could do to make it work?
I'm using QGIS 3.10.12


texts:
I´ve a project in VBA but I hit a wall that I can't pass.
What I want to do:

I´ve a table (of parcels) with the fields: ID, Name, Area, Scale. For each Parcel there is a different scale!

I input, in a text box, the name of a parcel. I want to "get" the corresponding scale.

Does anybody know how I can do that?


texts: Updating all attribute data in Feature Layer using ArcObjects

texts:
I have a layer with over 250,000+ features, and I want to update all the attribute data using ArcObjects and C#. I have written the code below and it takes around 3 hours to perform the update. Can you suggest whether the code below is efficient enough and any suggestion regarding the delta tables please?

public void updateWithCursor(IWorkspace workspace, IFeatureClass featureClass, List<Update> updateList)
    {
        using (ComReleaser comReleaser = new ComReleaser())
        {
            IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)workspace;

            // Start an edit session
            workspaceEdit.StartEditing(false);

            // Start an edit operation
            workspaceEdit.StartEditOperation();
            // Use IFeatureClass.Search to create a search cursor.
            IFeatureCursor searchCursor = featureClass.Search(null, false);
            comReleaser.ManageLifetime(searchCursor);

            // Find the positions of the fields used to get and set values.
            int idIndex = featureClass.FindField("ID");
            int serialIndex = featureClass.FindField("SERIAL");
            int cityIndex = featureClass.FindField("CITY");
            int updateDateIndex = featureClass.FindField("UPDATEDATE");

            int counter = 0;
            IFeature feature = null;
            while ((feature = searchCursor.NextFeature()) != null)
            {
                string searchValue = feature.get_Value(idIndex).ToString();

                if (searchValue != "")
                {
                    Update upd = updateList.SingleOrDefault(x => x.id == searchValue);

                    if (upd != null)
                    {
                        feature.set_Value(serialIndex , upd.serialNo);
                        feature.set_Value(cityIndex , upd.city);
                        feature.set_Value(updateDateIndex, DateTime.Now);
                        feature.Store();
                    }
                    Console.WriteLine(counter++);
                }
            }

            workspaceEdit.StopEditOperation();

            // Stop the edit session
            workspaceEdit.StopEditing(true);
        }
    }


texts:
I am repeatedly getting values in the 1000s for a nearest neighbour analysis in QGIS 3.4.

I've tested a different dataset and getting the expected result (i.e. between 0 and 2.5). But from the dataset I'm using I'm getting huge numbers:

Observed mean distance: 9.579807163911296
Expected mean distance: 0.0002003651354400665
Nearest neighbour index: 47811.747003144774
Number of points: 48
Z-Score: 633690.2583200812

A screenshot of the vector layer's points and the applied filter:

[1]: https://i.stack.imgur.com/etZkd.png [1]
  [2]: https://i.stack.imgur.com/PEVJA.png [2]

The dataset can be accessed this OneDrive link: https://isbdragons-my.sharepoint.com/:u:/g/personal/ofidler_isb_bj_edu_cn/EXC4adHTCl5KlSP4UQNblYwBFb-7rO9-VLK367Uevgf_pg?e=LShz0z

My suspicion is QGIS isn't bounding the extent of the analysis properly. I've tried 'selected features only', which generated almost the same numbers. When I 'zoom to layer' it's the same extents as shown in the screenshot above.


texts: How to create DEM from contours without ArcGIS Spatial Analyst licence?

texts:
I'm trying to create a 3d DEM/DTM from a flat contour shapefile which has height data attached to each contour, but I'm not sure how to do this?

I've found this post, on creating a dem from contours using the Topo to Raster tool, (using the ANUDEM technology,) which looks great: [see: creating a dem from contours] but I don't have a Spatial Analyst licence, so wondering if there is another way to do this, using a standard ArcView licence (version 9.2)?

I'm currently investigating QGIS as an alternative, [via this process: How to generate a DEM from a contour Shapefile?] though I've only just downloaded QGIS, so while I'm a GIS newbie I am SLIGHTLY more familiar with Arc!

Any alternatives, help, or tips would be much appreciated!


texts: Http error: 404 after Install JAI and JAI Image on GeoServer with export _POSIX2_VERSION=199209

texts:
I was installing JAI and JAI Image, when I follow the documentation from Geoserver, the JAI work ok, without problems, but when I install JAI Image with the commands, from jdk8 folder (location) as the documentation from geoserver:

$ export _POSIX2_VERSION=199209
$ sh jai_imageio-1_1-lib-linux-i586-jdk.bin
accept license_

the install is correct, but when i start Geoserver, it show the error message:

  http error: 404  problem accessing /geoserver/.  reason: not found

I can't access to Geoserver, i have to re install Java (Rollback), then the Geoserver working correctly.


texts: How to make a spatially map using interpolation method?

texts:
I have a .csv file which has PH value of groundwater for the couple of locations in state of WEST BENGAL (India). First I uploaded the .csv file data using "Add Delimited Text Layer" and make a Continuous distribution raster layer which have five classes and look like this:

Now I want to distribute these PH values of five classes in the Polygon shapefile of West Bengal state.  But when I import Polygon shapefile of West Bengal then the coordinates do not match with the raster field.

What should I do?


texts: Merging DEMs or other raster datasets and clipping Raster dataset to vector line?

texts:
I have two DEMs that I am going create a hillshade from, and ultimately want one hillshade clipped to a county boundary. I am running into a couple problems, though.

When I try to merge the two hillshades created from the DEMs, ArcMap crashes. When I try to merge the two DEMs, ArcMap crashes. I am assuming that if I managed to get the DEMs merged, the hillshade created would cause ArcMap to crash, too.

I am using the PLTS Merge Raster Datasets tool, because for some reason the Spatial Analyst tool Raster Calculator does not appear as selectable.

Once I am able to do this and get a hillshade of the two DEMs, is there a specific function for using the county line (which is a vector fc) to clip the raster data set to?


texts: Adjust label/symbol size in print composer automatically when choosing different scale for atlas

texts:
I am having trouble with my labels/symbols in the print composer and didn't find a solution here so far.

I have a project with a certain area from which I have to create several layouts in the print composer. Example below.

Map with overview in 1:250
Map with overview in 1:150 with several map sheets which are 1:50 or 1:20

When I create the map sheets with the atlas function in 1:50 and 1:20 it's necessary to adjust the label and symbol size for the layout in the print composer to get a good symbol/label size.

First image 1:20 and second 1:50

Is it possible to create a label or symbol rule which specifies the exact size of the symbol/label in the print composer when choosing a map sheet with 1:20 and 1:50? At the moment I manually change the size of the symbol/label for the print composer which is annoying.

Summarized:

map sheets created which are in 1:20 and 1:50 and show different things
the label/symbol size should be changed automatically when choosing 1:20 or 1:50 in the print composer via the atlas

Would this be possible?


texts:
I am trying to perform the kriging prediction of my data.
I follow a pretty standard process.
Firstly, I prep a grid to which the kriging predictions are supposed to be inserted.
I prep data that I want to interpolate.
I make sure that CRS is the same for both of the objects. Then as a next step, the variogram is selected. Lastly, we interpolate the values using kriging.
Yet, for some reason, the returned kriging object contains all predictions that are NAN.
I tried different CRS systems and I am not able to come up with a solution to this problem.
Therefore, I would like to ask for any suggestions.
Example Here:
library(tidyverse)
library(sf)
library(raster)
library(stars)
library(gstat)
library(ggmap)

# get data. using your example, we'll take spatial data from the raster package
Regions <- getData("GADM", country = "CZ", level = 1)

Regions <- Regions[Regions$NAME_1 == 'Prague', ]

Regions <-
  Regions %>%
  st_as_sf()

grid_spacing <- 0.005

# Create a grid for the border --------------------------------------------

polygony <-
  st_make_grid(Regions,
               square = T,
               cellsize = c(grid_spacing, grid_spacing)) %>%
  st_sf()

grid <- st_intersection(polygony, Regions)

plot(grid)

grid2 <- as_Spatial(grid$geometry)
plot(grid2)

coordinates(grid2)

spdf <- grid2
grd <- makegrid(spdf, n = 3000, pretty = F)

colnames(grd) <- c('x','y')

grd %>% plot()

grd_pts <- SpatialPoints(coords = grd,
                         proj4string=CRS(proj4string(spdf)))

# find all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]
grd_pts_in %>% plot()

# transform grd_pts_in back into a data frame
gdf <- as.data.frame(coordinates(grd_pts_in))

gdf %>%
  mutate(value = rnorm(nrow(gdf))) %>%
  ggplot(aes(x, y, fill =  value)) +
  geom_tile()

grid2 %>% plot()
gdf %>% plot()

system <- "+proj=utm +zone=19 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
system2 <- "+proj=longlat +datum=WGS84 +no_defs "

# Data to be interpolated -------------------------------------------------

df <-
  tibble(
long = sample(coordinates(grid2)[,1], 500),
lat = sample(coordinates(grid2)[,2], 500),
val = rnorm(500, 5000, 10)
)

df %>%
  ggplot(aes(long, lat, color = val)) +
  geom_point()

# Kriging Predictions Itself ----------------------------------------------

resid2 <- df %>%
  st_as_sf(coords = c("long", "lat"), crs = system2)

cutoff = 3
lzn.vgm <- variogram(val ~ 1, data = resid2, cutoff = cutoff)
plot(lzn.vgm)

nugget <- 100 # is the value of the semivariance at zero distance.
psill <- max(lzn.vgm$gamma) - nugget # is the difference between the sill and the nugget
range <- 0.2 # distance where the model first begins to flattens out.

lzn.fit.Ste <- fit.variogram(lzn.vgm, model = vgm(
  model = "Ste",
  nugget = nugget,
  sill = psill,
  range = range))

plot(lzn.vgm, lzn.fit.Ste)

coordinates(grid2)

POLY2 <- as(grid2, "Spatial")
RES <- as(resid2, "Spatial")

crs(RES) = CRS(system2)
crs(grid2)  = CRS(system2)

st_crs(grid2) == st_crs(RES)

OK =
  krige(val ~ 1, RES, grid2,
        model = lzn.fit.Ste)

OK %>%
  as.data.frame()


texts: Counting points in a scratch layer polygon using QGIS

texts:
I'm using QGIS 3 and I am trying to count the points in a scratch layer polygon. When I use the count points in a polygon tool it says there are no points in my polygon which is untrue.
I figured the polygon being a scratch layer could be the issue. I get an error message when trying to create the polygon in a shapefile. So i'm not sure how to fix that either.


texts: How to make rectangle selection on QgsMapCanvas?

texts:
I have this code for create and open raster map:

if(QFile::exists(rasterFileName))
{
    QFileInfo myRasterFileInfo(rasterFileName);
    m_pRasterLayer = QSharedPointer<QgsRasterLayer>(
        new QgsRasterLayer(myRasterFileInfo.filePath(),
        myRasterFileInfo.completeBaseName()));

    if (m_pRasterLayer->isValid())
    {
        QgsMapLayerRegistry::instance()->addMapLayer(
            m_pRasterLayer.data(), true);

        QList<QgsMapCanvasLayer> layerList;
        layerList.append(QgsMapCanvasLayer(m_pRasterLayer.data(), true));

        m_pMapCanvas->setLayerSet(layerList);

        /* ... */
    }
}

How can I make (write code for) rectangle selection by mouse on my map?


texts: I can't find the New label tools in QGIS developed by sourcepole. I'm using QGIS 1.7.4. Is it a QGIS tool or a plugin?


texts: I have an *.mxd file and I don't want to save any edits to it after each session. Is there a way in 'arcpy' or python in general to restrict edit of *.mxd file? Any suggestions would be appreciative.


texts:
My python console won't work at all. I can't open it.
My QGIS version is 2.4. OS is win7 64.
I've reinstalled QGIS a few times but I always get the same result.
Strange. Something went wrong with the source code I guess?
Here's the error :

Traceback (most recent call last):   File "", line 2, in    File "C:/OSGEO4~1/apps/qgis/./python\console\console.py", line 43, in
 show_console
     _console = PythonConsole( parent )   File "C:/OSGEO4~1/apps/qgis/./python\console\console.py", line 75, in
 __init__
     self.console = PythonConsoleWidget(self)   File "C:/OSGEO4~1/apps/qgis/./python\console\console.py", line 101, in
 __init__
     self.shellOut = ShellOutputScintilla(self)   File "C:/OSGEO4~1/apps/qgis/./python\console\console_output.py", line 103,
 in __init__
     self.insertInitText()   File "C:/OSGEO4~1/apps/qgis/./python\console\console_output.py", line 145,
 in insertInitText
     socket.gethostname()) UnicodeDecodeError: 'ascii' codec can't decode byte 0xe9 in position 9: ordinal not in range(128) Version de
 Python :
 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)] Version de QGIS :
 2.4.0-Chugiak 'Chugiak', 8fdd08a Chemin vers Python : ['C:/OSGEO4~1/apps/qgis/./python/plugins\\processing',


texts: Finding the coordinate system info from a LAS/LiDAR file with Python?

texts:
I found several examples of how one can query and manually set a coordinate system on a tif/raster file (GDAL/Python: How do I get coordinate system name from SpatialReference? or Is it possible to get the EPSG value from an OSR SpatialReference class using the OGR Python API? for example). But if I have a LAS file with Lidar point data and want to read its coordinate system programmatically with python, how could I go about doing that? I tried the two scripts below with no luck:

# first:

import gdal, ogr, os, osr
from laspy.file import File
driver = ogr.GetDriverByName('LOSLAS')  ## maybe this driver would read LAS info, but hasn't worked
dataset = driver.Open(r'file.las')  ## this could be the problem
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef = geom.GetSpatialReference()
print spatialRef, geom

# second:

from osgeo import gdal,osr
dataset = gdal.Open(r'file.las')  ## again this command doesn't work
spatialRef = osr.SpatialReference(wkt = dataset.GetProjection())  ## get projection function also throws an error
print 'projection = ', spatialRef
if spatialRef.IsProjected:
    print spatialRef.GetAttrValue('projcs')
print spatialRef.GetAttrValue('geogcs')

Any thoughts on what I am doing wrong or what other libraries I could try? Again I want to be able to read the spatial referencing of a LAS file using python, without having to convert it to a raster first or use LAStools functionalities. Maybe there's a way to query this within the header?


texts:
With R I want to make an animated choropleth map with the rMaps package  and a custom map. I tried to follow the example of Diego Valle described here and to use the custom map info for Datamaps.

However, the map does not display. But, the legend with the proper fill key appears with the year animation bar. The topoJson file seems to be fine.

library(rCharts);library(rMaps);library(rgdal)
map = readOGR("https://raw.githubusercontent.com/gupierre/GHA_vol/master/Data/GHA_adm/GHA_adm1.topo.json", layer="admin_1",stringsAsFactors=F)
#Data:

SDLOG<-read.table("https://raw.githubusercontent.com/gupierre/GHA_vol/master/Data/sdlog.csv", header=T,stringsAsFactors=F, sep =",")
head(SDLOG)
 year          name      sdlog  type
>1 2004 Greater Accra 0.20030687 urban
>2 2005 Greater Accra 0.15322712 urban
>3 2006 Greater Accra 0.10782677 urban

#Map:

library(rCharts);library(rMaps)
d1 <- ichoropleth(sdlog ~ name, data = SDLOG[SDLOG$type=="rural",], ncuts = 7, pal = 'YlOrRd', animate = 'year',  map = 'admin_1')
d1$set(
  geographyConfig = list(
    dataUrl = "https://raw.githubusercontent.com/gupierre/GHA_vol/master/Data/GHA_adm/GHA_adm1.topo.json"
  ),
 scope = 'admin_1',
  setProjection = '#! function( element, options ) {
  var projection, path;
  projection = d3.geo.mercator()
  .center([8, -1]).scale(element.offsetWidth)
  .translate([element.offsetWidth / 2, element.offsetHeight / 2]);

  path = d3.geo.path().projection( projection );
  return {path: path, projection: projection};
  } !#'
)
d1

Note: this question was asked on stackoverflow without success. I don't know the policy on GIS for posting questions from SO. Please advise in the comments if I am wrong in simply copy pasting my initial question.


texts: How to add an attribute with QGIS 2.8.1 Python console?

texts:
i would like to add an attribute with the python console:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
layer = iface.addVectorLayer("C:/donnees_SIG/05.shp","nou","ogr")
layer = iface.activeLayer()
layer.startEditing()
layer.addAttribute(QgsField("STATUT", QVariant.Int))

But it doesn't work...
When i write the last line, the python console returns 'False'.

I'm using QGIS 2.8.1


texts: Adding alpha band to a TIFF image drastically increases its filesize?

texts:
I'm using gdalwarp -dstalpha src.tif dst.tif to add alpha channel to an 8-bit unsigned image with 3 bands. The initial size of the image is 300 mb, but adding an alpha band increases it to 1980 mb, which is far more than expected 1/4 increase.

What is the matter?


texts: Getting an error when trying to set extent for raster() in R

texts:
I'm working through Robert's answer to this question on how to plot irregular raster data: How to make RASTER from irregular point data without interpolation
no2tc <- ncvar_get(nc, "DETAILED_RESULTS/nitrogendioxide_total_column")
lat <- ncvar_get(nc, "PRODUCT/latitude")
lon <- ncvar_get(nc, "PRODUCT/longitude")

no2vec <- as.vector(no2tc)
latvec <- as.vector(lat)
lonvec <- as.vector(lon)

df <- data.frame(lonvec,latvec)
colnames(df) <- c('X', 'Y')
e <- extent(df)

My data very similar in format to the fake data created in his question and I got the following error:
Error in extent(df[, 2:3]) : c("x", "y") %in% names(x) are not all TRUE

Here is a sample of my data
     Z         X          Y
1   NA -84.17967 -118.36866
2   NA -84.22738 -117.62881
3   NA -84.27346 -116.88749
4   NA -84.31795 -116.14489
5   NA -84.36089 -115.40117
6   NA -84.40232 -114.65651
7   NA -84.44226 -113.91109
8   NA -84.48075 -113.16512
9   NA -84.51782 -112.41875

> df[,2:3]
            X             Y
1   -84.17967 -118.36865997
2   -84.22738 -117.62880707
3   -84.27346 -116.88748932
4   -84.31795 -116.14488983
5   -84.36089 -115.40116882

The data has some NA values but obviously, I would want to ignore these.
Can someone tell me why I might be getting this error?


texts:
I am making a legend in ArcGIS 10, and have a few layers with a large number of classes. When I try to display the legend with two columns for that layer, the second layer always lines up with the heading/layer name. This doesn't look proper, and I want the first item in each column to line up with each other. I have used a few work-arounds, like splitting the legend, or converting to graphics and manually adjusting.
Is there a place in the legend properties where this can be done?


texts: Changing SQL syntax in QGIS Virtual Layers

texts:
I am experiencing difficulties with SQL syntax in a Virtual Layer. I can do the simplest query form such as:

SELECT *
FROM mytable

and this displays on my canvas.

However, if I try to name fields in my table such as:

SELECT field1, field2
FROM mytable

this doesn't create an error but does not display on my canvas.

Can anyone offer suggestions, please? I have tried it on several layers.


texts: Function fires despite conditional query in leaflet-map

texts:
In a map built with Leaflet I created a button to export and download data as geoJSON drawn before by the user with Leaflet.draw.

To prevent an empty download, I am querying the existence of data.

In the leaflet drawing tool, the created data is stored in a layergroup named 'layerGroup'.

After loading the map for the first time, by clicking on the button I am getting the alert "There is nothing to download!", the console returns "Layers in featureGroup: 0" and no data was downloaded. So far so good.

Within the drawing tool I am able to add and remove data later again. After adding and then removing all data, there should be no layer in layerGroup and therefore I should experience the same acting like before with no download. In the console I get the message "Layers in featureGroup: 0", the alert "There is nothing to download!" is showing up, but, I don't know why, the download still starts.

Do you have any ideas why my code is acting like that?

That's my current attempt:

//Export Button
document.getElementById('export').onclick = function(e) {
        //get layer number in var layer_number
        var layer_number = featureGroup.getLayers().filter( function(l) {return l instanceof L.Marker} ).length;
        console.log("Layers in featureGroup: ", layer_number);
        if (layer_number > 0){

            // Extract GeoJson from featureGroup
            var data = featureGroup.toGeoJSON();

            // Stringify the GeoJson
            var convertedData = 'text/json;charset=utf-8,' + encodeURIComponent(JSON.stringify(data));

            // Create export
            document.getElementById('export').setAttribute('href', 'data:' + convertedData);
            document.getElementById('export').setAttribute('download',''+timestamp +'.geojson');

            //show confirmation popup
            $('.confirmation').show();
            //hide popup on button click
            $('#agree_button').click(function(){
                $('.confirmation').hide();
            });
        }
        //meassage if no feature in featureGroup
        else {
            alert("There is nothing to download!");
        }
    }


texts: GeoServer and TIGER roads layer overlay problem

texts:
Been playing with OpenLayers, GeoServer and the TIGER roads shapefile but cannot get it to line up correctly over New York.  I have looked online and followed the instructions regarding setting the projections and setting Google Maps to sphericalmercator = true, however it still doesn't line up correctly.

Any ideas?

        var map;
        var untiled;
        var tiled;
        var pureCoverage = false;
        // pink tile avoidance
        OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;
        // make OL compute scale according to WMS spec
        OpenLayers.DOTS_PER_INCH = 25.4 / 0.28;

        function init(){
            format = 'image/png';

            var bounds = new OpenLayers.Bounds(
                -20037508, -20037508,
                20037508, 20037508.34
            );

            var options = {
                projection: new OpenLayers.Projection("EPSG:900913"),
                displayProjection: new OpenLayers.Projection("EPSG:4326"),
                maxExtent: bounds,
                //maxResolution: 0.00075764453125,
                units: 'm'
            };

           map = new OpenLayers.Map('map', options);

           var GoogleMap = new OpenLayers.Layer.Google(
                "Google Satellite", {
                type: G_SATELLITE_MAP,
                maxZoomLevel: 30,
                "sphericalMercator": true
                }
           );

            // setup tiled layer
            wmsLayer = new OpenLayers.Layer.WMS(
                "tiger_roads - Tiled", "http://localhost:8080/geoserver/wms",
                {
                    layers: 'tiger:tiger_roads',
                    styles: '',
                    format: format,
                    tiled: 'true',
                    transparent:'true',
                    isBaseLayer: false
                }
            );

            map.addControl(new OpenLayers.Control.MouseDefaults());
            map.addControl(new OpenLayers.Control.Navigation());
            map.addControl(new OpenLayers.Control.ZoomBox({alwaysZoom:true}));
            // Coordinate display at bottom of map
            map.addControl(new OpenLayers.Control.MousePosition());
            map.addLayers([GoogleMap, wmsLayer]);

           var point = new OpenLayers.LonLat(-73.96543, 40.78885);
           // Need to convert zoom point to mercator too
           point.transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());
           map.setCenter(point, 12);

      }


texts: Project one raster into another raster's crs in gdal

texts:
I've two rasters and would like to project one into coordinate system of another. One raster has an extent in meters and one in lat/long. I would like use gdal to take the raster which has an extent in lat/long and to project it into crs of another raster.

The output of gdalinfo of the first raster is:

PROJCS["NZGD_2000_New_Zealand_Transverse_Mercator",
    GEOGCS["GCS_NZGD_2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS_1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1323824.000000000000000,5190684.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 1323824.000, 5190684.000) (169d35'25.39"E, 43d23' 9.18"S)
Lower Left  ( 1323824.000, 5004684.000) (169d29'32.65"E, 45d 3'26.70"S)
Upper Right ( 1422924.000, 5190684.000) (170d48'45.84"E, 43d24'56.95"S)
Lower Right ( 1422924.000, 5004684.000) (170d44'59.23"E, 45d 5'20.90"S)
Center      ( 1373374.000, 5097684.000) (170d 9'43.02"E, 44d14'19.75"S)

And the output of the gdalinfo for the second raster is:
GEOGCS["NZGD49",
    DATUM["New_Zealand_Geodetic_Datum_1949",
        SPHEROID["International 1924",6378388,297.000000000005,
            AUTHORITY["EPSG","7022"]],
        TOWGS84[59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993],
        AUTHORITY["EPSG","6272"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4272"]]
Origin = (166.474999999999990,-34.424999999999997)
Pixel Size = (0.049382716049383,-0.049423076923077)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 166.4750000, -34.4250000) (166d28'30.00"E, 34d25'30.00"S)
Lower Left  ( 166.4750000, -47.2750000) (166d28'30.00"E, 47d16'30.00"S)
Upper Right ( 178.4750000, -34.4250000) (178d28'30.00"E, 34d25'30.00"S)
Lower Right ( 178.4750000, -47.2750000) (178d28'30.00"E, 47d16'30.00"S)
Center      ( 172.4750000, -40.8500000) (172d28'30.00"E, 40d51' 0.00"S)

Basically I need to repoject second raster into crs of first and have extent in meters.


texts: Changing Mosaic Dataset Properties in Map Document using ArcPy?

texts:
Does anyone know if it is possible in ArcGIS 10.3 for Desktop to access and change mosaic dataset properties within a map document using ArcPy?

I am trying to write text into the symbology description in the boundary layer of a mosaic dataset (not the description field on the layer properties general tab of the mosaic dataset).


texts: How to deactivate Vector Selection In OpenLayers wms features

texts:
This is my sample code and work fine :), But I have a little problem,

I'm using multiple layers (WMS), when I zoom  the layer will be shown or hidden.
In my example bellow when OpenLayers level zoom is greater than 3 I'm showing layer2 and hiding layer1. And also I create a control to active highlight ability.
this example code for

        control = new OpenLayers.Control.GetFeature({
        protocol: OpenLayers.Protocol.WFS.fromWMSLayer(layer2),
         box: true,
         multipleKey: "shiftKey",
         toggleKey: "ctrlKey"
        });

        control.events.register("featureselected", this, function(e) {
              select.addFeatures([e.feature]);
        });
        control.events.register("featureunselected", this, function(e) {
              select.removeFeatures([e.feature]);
         });

        map.addControl(control);
        control.activate();

when OpenLayers level zoom less than 3 I'm showing back layer1 and hiding layer2. and I will deactivate highlighting control. and clearing all features selected.

I have trouble to do it, and bellow example code :

    // Layer 1
    layer1 = new OpenLayers.Layer.WMS(
        "Map A", serverUrl,
        {
            LAYERS: 'basic:layerA',
            STYLES: '',
            format: format
        },
        {
            isBaseLayer: true,
            visibility: true
        }
    );

    // Layer2
    layer2 = new OpenLayers.Layer.WMS(
        "Map B", serverUrl,
        {
            LAYERS: 'basic:layerB',
            STYLES: '',
            format: format,
            transparent: true,
            tiled: true,
        },
        {
            isBaseLayer: false,
            visibility: false,
            unsupportedBrowsers:[]
        }
    );

    select = new OpenLayers.Layer.Vector("Selection", {styleMap:
        new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])
    });

    map.addLayers([layer1,layer2,select]);

//Zoomend event
    map.events.on({"zoomend": function (e) {

        console.log(this.getZoom());

    //I have trouble in here
        if (this.getZoom() > 3) {
            layer2.setVisibility(true);
            layer1.setVisibility(false);

            control = new OpenLayers.Control.GetFeature({
            protocol: OpenLayers.Protocol.WFS.fromWMSLayer(layer2),
             box: true,
             multipleKey: "shiftKey",
             toggleKey: "ctrlKey"
            });

            control.events.register("featureselected", this, function(e) {
                select.addFeatures([e.feature]);
            });
            control.events.register("featureunselected", this, function(e) {
                select.removeFeatures([e.feature]);
            });

            map.addControl(control);
            control.activate();

        }
        else{
            // In here, I will deactive controll above (hightligthing controll)

            layer2.setVisibility(false);
            layer1.setVisibility(true);
        }
    }
    });

And should must I do ?


texts:
I want to calculate the area for snow cover but it doesn’t work out for me.

var collection = ee.ImageCollection('COPERNICUS/S2')
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10))
  .filterDate('2020-04-01' ,'2020-04-27')
  .filterBounds(polygon) ;

print(collection) ;

  //Let's centre the map view over our ROI
Map.centerObject(polygon, 13);
var medianpixels = collection.median(); // This finds the median value of all the pixels which meet the criteria.

var medianpixelsclipped = medianpixels.clip(polygon).divide(10000);
// Now visualise the mosaic as a natural colour image.
Map.addLayer(medianpixelsclipped, {bands: ['B4', 'B3', 'B2'], min: 0, max: 1, gamma: 1.5}, 'Sentinel_2 mosaic',false);

// User specified parameters

// setup Visualization
var viz = {min:-0.50, max:0.75, palette:['red', 'green', 'orange', 'white']};

// Data
// calculate ndsi
var ndsi = medianpixelsclipped.normalizedDifference(['B3', 'B11']);
var scale=30;
Map.addLayer(ndsi, {min:0.4, max:1, palette:['EB984E ','FDFEFE  ']}, 'ndsi');
var reducer = ndsi.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: polygon,
scale: 30,
maxPixels: 1E13
});
// km square
var area = ee.Number(reducer.get('NDSI')).multiply(scale).multiply(scale).divide(1000000);
print('area of ndsi ', area.getInfo() + ' km2');


texts: QGIS 2.2 Print Composer does not display map with some datasets

texts:
I can't seem to get the Print Composer of QGIS 2.2 to work. When I click on "Add new map", the map square is created, but I can't seem to see my data in it, whatever the scale I use or even clicking on "set to map canvas extent". When I add a Legend, the contents are correct, all the layers I have rendered on the main QGIS window are listed there.

Can anyone help me with this, please? Thanks!


texts: I am working in a project that involves data from different files types, I managed to load TIFF and .Gz file (ordinance survey mastermap topography) to correct location. However I have a DXF layer that also needs to be overlayed but when I import it, the location is nowhere near the other two layers. As far as I can tell all my CRS settings are correct.  Does anyone have any ideas what I can do to fix this issue?


texts: Displaying annotation in older versions of ArcGIS Desktop?

texts:
I created an annotation layer in a trial of 10.5 that will not open in 10.0. I don't see why this should be as the annotation resides in my fgdb and not a map document.

When I try to open it in 10.0 I get the following error:

"Item does not have a definition. The XML being loaded could not be parsed."

The associated layer opens fine.


texts:
Regarding my pyqgis plugin 2.2 , I've implemented signal on a vector layer event  successfully. For instance, whenever i create a new feature, i can read a new basic message appeared via QGIS logging system : " feature added".

But for some  events, besides this basic message, i would like read the name of the layer and/or , the id of the feature/ the attribut name / the geometry type (line,surf,punct).

How to bring more info details for the events?

I paste below my piece of code :

def logLayerModified(self, onlyGeometry = None ):
    QgsMessageLog.logMessage( "layer modified" )
    QApplication.beep()

def logAttributeAdded(self,QgsFeatureId):
    QgsMessageLog.logMessage( "attribut added" )
    QApplication.beep()

def logAttributeDeleted(self,QgsFeatureId):
    QgsMessageLog.logMessage( "attribut deleted" )
    QApplication.beep()

def logFeatureAdded(self,QgsFeatureId):
    QgsMessageLog.logMessage("feature added, id = " + str(QgsFeatureId))
    QApplication.beep()

def logFeatureDeleted(self,QgsFeatureId):
    QgsMessageLog.logMessage("feature deleted, id = " + str(QgsFeatureId))
    QApplication.beep()

def logGeometryChanged(self,fid):
    QgsMessageLog.logMessage("Geometry changed, id = " + str(fid))
    QApplication.beep()

def logAttributeValueChanged(self, QgsFeatureId fid, int idx, const QVariant &):
    QgsMessageLog.logMessage("Attribut value changed", id = " + str(QgsFeatureId))
    QApplication.beep()

def logBeforeLayerCommitChange(self):
    QgsMessageLog.logMessage( "commit change" )
    QApplication.beep()

layer.layerModified.connect(self.logLayerModified)
layer.attributeAdded.connect(self.logAttributeAdded)
layer.attributeDeleted.connect(self.logAttributeDeleted)
layer.featureAdded.connect(self.logFeatureAdded)
layer.featureDeleted.connect(self.logFeatureDeleted)
layer.featureDeleted.connect(self.logGeometryChanged)
layer.beforeCommitChanges.connect(self.logBeforeLayerCommitChange)
layer.AttributeValueChanged.connect(self.logAttributeValueChanged)