texts: Converting Google Maps lat long to EPSG:6706

texts:
I'm sending a WMS GetFeatureInfo() request, but there is something wrong because the returned info does not appropriate a clicked point.

I suspect that the error is in the BBOX parameter of the request because I'm not sure that bounds are in the correct EPSG.

How can I convert Google Maps Bounds to the EPSG:6706?

This is what I tried to do:

google.maps.event.addListener(map, "dblclick", function (e) {
  if ($('#wms').is(":checked")) {
  var scale = Math.pow(2, map.getZoom());
  var proj = map.getProjection();
  var bounds = map.getBounds();
  var nw = proj.fromLatLngToPoint(
      new google.maps.LatLng(
      bounds.getNorthEast().lat(),
      bounds.getSouthWest().lng()
  ));
  var so = proj.fromLatLngToPoint(
      new google.maps.LatLng(
         bounds.getNorthEast().lng(),
         bounds.getSouthWest().lat()
  ));
  var point = proj.fromLatLngToPoint(e.latLng);
  var i = Math.floor((point.x - nw.x) * scale);
  var j = Math.floor((point.y - nw.y) * scale);
  var w = $('#map').width();
  var h = $('#map').height();
  i = Math.round(i*256/w);
  j = Math.round(j*256/h);
  var url = "https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?";
  url += "SERVICE=WMS";    //WMS service
  url += "&VERSION=1.3.0";  //WMS version
  url += "&REQUEST=GetFeatureInfo"; //WMS operation
  url += "&FORMAT=image/png" ; //WMS format
  url += "&TRANSPARENT=TRUE";
  url += "&QUERY_LAYERS=CP.CadastralParcel"; //,CP.CadastralZoning";
  url += "&LAYERS=CP.CadastralParcel"; //,CP.CadastralZoning";
  url += "&TILED=true";
  url += "&INFO_FORMAT=text/html";
  url += "&I="+i;
  url += "&J="+j;
  url += "&WIDTH=256";         //tile size in google
  url += "&HEIGHT=256";
  url += "&CRS=EPSG:6706";     //set WGS84
  url += "&BBOX=" + bbox;      // set bounding box
  window.open(url);
}
});

This is an example of the WMS call: it works but the returned info are not correct for the point clicked.

https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=TRUE&QUERY_LAYERS=CP.CadastralParcel&LAYERS=CP.CadastralParcel&TILED=true&INFO_FORMAT=text/html&I=65&J=88&WIDTH=256&HEIGHT=256&CRS=EPSG:6706&BBOX=10.510482788085938,44.35429627478576,10.51116943359375,44.35478724525106


texts: QGIS 3.0 Girona using profile tool, get python error

texts:
I have a raster layer (LIDAR data) selected. On releasing the mouse button, and before any rubber band starts, I get this error:
Message log:

2018-03-09T16:19:08 WARNING Traceback (most recent call last):
              File "C:/Users/John/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\profiletool\tools\ptmaptool.py", line 78, in moved
                self.profiletool.rubberband.addPoint(QgsPoint(self.pointstoDraw[i][0],self.pointstoDraw[i][1]))
            TypeError: QgsRubberBand.addPoint(): argument 1 has unexpected type 'QgsPoint'

Can you help, please?


texts: Generating urban footprint polygons from building clusters in QGIS?

texts:
I'm not sure what the correct term is, but it's a common cartographic design pattern.

I want to show the outline of built-up areas as a filled polygon beneath road and buildings layers.

In OpenStreetMap, such areas are sometimes available as landuse=residential. I suspect these will all have been traced manually.

Here's an example, this is Gibraltar. I digitzed this outline manually and clipped it to the shoreline. The "urban footprint" consists of the purple/pink areas (this was extended out to sea, then clipped to the land outline from OpenStreetMapData)

  Is there a tool available in QGIS (or other FOSS GIS tools) to generate these automatically from a set of buildings?

Two approaches I've tried already...

I tried the Concave hull plugin but only got a single area.
I also tried a Delauney Triangulation on nodes extracted from the building outlines. I suspect if I can discard large polygons or those with very acute angles between vertices, that might work...


texts:
I'm trying to take a max annual NDVI product derived from AVHRR and relate it to a specific area that I'm examining.This requires me to identify the region of interest based on Landsat scenes, and relate it to AVHRR's raster resolution. I know what area I want to look at in the Landsat scenes, but am having trouble translating it to area in AVHRR, which is currently projected in North Pole Lambert Azimuth Equal Area.

My landsat layer is currently projected in UTM Zone 13N.

I have tried switching both projections to something compatible, with each to USA Contigous Lambert Conformal Conic, but there still seems to be a fundamental projection issue occurring. The region of interest is in NM, and it doesn't appear the AVHRR scene is lining up as shown:

The basemap is projecting in WGS84, and it appears as though the reprojected landsat scene conforms with the baseman, whereas the AVHRR scene does not.


texts: Managing coordinates and projections for zonal statistics on global data

texts:
I am working with a global raster that is in the World Mollweide (EPSG: 54009) projection. I also have a set of polygons from different parts of the globe that are shapefiles in the WGS 84 (EPSG: 4326) CRS.
I am well aware that the best approach for summing the raster scores within each polygon would be to project the shapefiles and then conduct zonal statistics with both layers in Mollweide (since projections between rasters involve resampling of the grid cell values). However, the zonal statistics will be part of a larger geoprocessing setup where the conversion of the polygons/shapefiles to Mollweide may not be as feasible as having the raster data in WGS 84 right from the start.
I was therefore wondering if there would be an additional way to conduct zonal statistics with both layers in WGS 84 while minimizing the errors that result from converting the raster files from Mollweide to WGS 84?
I have conducted a few analyses on how the results would differ if I were to use both the raster + shapefiles in Mollweide vs. having them both in WGS 84. What I found is that, relative to Mollweide, WGS 84 tends to overestimate the sums of scores within each polygon for locations that are further away from the equator while it underestimates the sums for locations that are close to the Equator.
If is it helpful, I am happy to provide some sample data and more details on how I compared the sums with both data inputs.


texts: Extracting projection details for TAB files from multiple folders using arcpy's describe spatial reference?

texts:
I am trying to extract projection details of vector files (it might be a tab or shp file) in multiple subfolders using arcpy's describe.SpatialReference. If i run this code, i am getting this error. I've just changed the code given here from list feature_class to list files. Since i am very much new to arcpy, don't know what went wrong

"AttributeError: DescribeData: Method spatialReference does not exist "

import arcpy
import os

arcpy.env.workspace = "D:/QAQC"
feature_classes = arcpy.ListFiles()

for csv_file in arcpy.ListFiles("*.tab"):
spatial_ref = arcpy.Describe(fc).spatialReference
  if spatial_ref.name == "Unknown":
     print("{0} has an unknown spatial reference".format(fc))
  else:
     print("{0} : {1}".format(csv_file, spatial_ref.name))

Edit: I've run the same code for shp and it's working perfectly.


texts: Calibrate routes not lining up with calibration points

texts: I am currently building a model that calibrates Z and M values from Cross Sections for a stream. Z being the 100yr flood and M being the station number.  What's happening is the vertex on the stream where the cross section intersects does not contain the 100yr flood and station number of that cross section.  For the picture below, the red is the true value of the cross section and the blue is where the cross section value should be.


texts: Recalculating area of single field in another projection?

texts:
I have polygon area. ArcGIS automatically calculate it's area for current coordinate system. I would like to create second field and calculate it's area in another coordinate system.

So in result I would like have two fields with two different coordinate system.

Is it possible to do in ArcGIS? Only for field, no changing projection of all data.


texts: How to get data from table and define the geometry on the same query?

texts:
My team is using DB2 database and I have to create and publish a map service in ArcGIS using one of the tables. The table has the XYZ column but no SRID defined in the table or database. The DB2 spatial extender is installed.

I tried to import using the new query layer but the geometry cannot be defined and I end up with an "Event table" layer which cannot be published.

I guess that there is a way to define the geometry via the query, but I am not sure how to do so. So, the question is: how can I retrieve data from table (X,Y,Z) and define the geometry on the same query?


texts: Difficulty with select by attributes using list fields in for loop

texts:
I am trying to add fields to a table and then populate the fields with text values (Yes/No). The text values are based on whether or not a species was found in a specific block. If the species was found in a block then the field > 0 and if no species was found in a block then field = 0.

I am having issues with my select by attribute where clause. Specifically, I am trying the use the field name text from my ListFields function in the sql statement (where_clause= field4 + " > 0").

I get the following error:

  ExecuteError: ERROR 000358: Invalid expression

import arcpy
from arcpy import env

arcpy.env.overwriteOutput = True

# Purpose is  to generate a table with presence/absence data for blocks falling outside of species ranges

blocks = r"C:\Users\Range\Workspace\scrap\FB_Ground_MLR.shp"

# layer for make feature

blyr   = r"C:\Users\Workspace\scrap\FB_Ground_MLR.lyr"
lyr = arcpy.MakeFeatureLayer_management(in_features = blocks, out_layer = blyr)

fields = arcpy.ListFields(blocks)
x = 0

for field in fields:
    field4 = field.name
    f4 = field4[0:8] + "X"
    x = x + 1
    if x >= 4:
    arcpy.AddField_management(blocks, f4, "TEXT")
    arcpy.SelectLayerByAttribute_management(lyr, selection_type="NEW_SELECTION", where_clause= field4 + " > 0")
    arcpy.CalculateField_management(blocks, f4, 'Yes')
    arcpy.SelectLayerByAttribute_management(lyr, selection_type="NEW_SELECTION", where_clause= field4 + " = 0")
    arcpy.CalculateField_management(blocks, f4, 'No')
else:
    print(f4)


texts: Using Loop function with arcpy.Exists()?

texts:
I have a code snippet using elif that checks if a file exists and continues to check through the letters of the alphabet (for example _A, _B, _C etc). All works fine but I think there must a better (or cleaner) way to write the code using a loop function instead of the way I have done.

#Create MXD file path
mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01.mxd"
mxdPath = wrFolder + "/M_" + WorkrequestRaw +"_01"

#check for latest version and create next letter
if  not arcpy.Exists(mxdOutput):
    mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01.mxd"
    pdfversion = ""
elif not arcpy.Exists(mxdPath +"_A.mxd"):
    mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01_A.mxd"
    pdfversion = "A"
elif not arcpy.Exists(mxdPath +"_B.mxd"):
    mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01_B.mxd"
    pdfversion = "B"
elif not arcpy.Exists(mxdPath +"_C.mxd"):
    mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01_C.mxd"
    pdfversion = "C"
elif not arcpy.Exists(mxdPath +"_D.mxd"):
    mxdOutput = wrFolder + "/M_" + WorkrequestRaw +"_01_D.mxd"
    pdfversion = "D"

I am a self taught beginner and trying to keep the code tidy. Welcome any comments on a better methodology?

Update poorly written question: What I am trying to do is if file name exists write new file called 'filenameA' and if 'filenameA' exists write 'filenameB' and if 'filenameB exists write 'filenameC' etc.


texts: Modify the Geoserver REST response href values to use domain name instead of IP address

texts:
Currently, the Geoserver REST API responds with objects that look like such:

{
    "layers": {
        "layer": [
            {
                "name": "layername1",
                "href": "http://localhost:8080/geoserver/rest/layers/layername1.json"
            },
            {
                "name": "layername2",
                "href": "http://localhost:8080/geoserver/rest/layers/layername2.json"
            }
        ]
    }
}

Can I configure geoserver to use a domain name for the href values like so:

{
    "layers": {
        "layer": [
            {
                "name": "layername1",
                "href": "http://example.domain.com/geoserver/rest/layers/layername1.json"
            },
            {
                "name": "layername2",
                "href": "http://example.domain.com/geoserver/rest/layers/layername2.json"
            }
        ]
    }
}


texts: Adding graphs, pie charts, bar diagrams in geoexplorer?

texts:
Does any one know how to add interactive graphs in Geoexplorer as shown in this video ?

Do we have to incorporate Google fusion tables ?

https://www.youtube.com/watch?v=d9ADUHGUXL0


texts: Installing packages into OSGeo4W64-QGIS-Python and not in Standard Python installation on Windows?

texts:
On my Windows 10 PC there are (at least) 2 Python 3.7 versions installed:

a) "Standard Python 3.7 64 bit with IDLE, which I downloaded from the site python.org and installed.

b) The Python 3.7 which comes bundled with the QGIS network-installer OSGeo4W which installed QGIS in the folder c:\OSGeo4W64

Problem:

I need to install some Python-Packages via pip in the Python-version, which is used by QGIS. I dont want to "mess up" my Standard installation of Python.


texts: Visualizing map from Spatialite DB on Mobile

texts:
I am working on an Android application for a school project.

I am converting shapefiles (examples) into SQLite databases using SpatiaLite and performing queries on those databases on the mobile with SpatiaLite-Android.

I am still looking for a way to visually display the map data on Android for user interaction. For example, I would like the user be able to zoom in/out and pan the map around. I also intend show the user her current location with respect to the map.

It would be great if my app can draw maps directly from the SQLite database with SpatiaLite. Can anyone offer tips on this? If this is not possible or not advisable, can anyone suggest alternatives?

NOTE: My app needs to function in offline mode.


texts:
I want to create a vector tile base map using my own shapefiles/postgis tables for use on a mobile primarily. By my understanding the process is divided into 3 main steps:

Creating/slicing the data into vector tiles.
Using a server to host the tiles.
Rendering the tiles.

I have looked into various methods to do this, The Mapzen Vector Tiles were the first thing I looked into. They come with a very good rendering option as well. I got it to work with osm data but I wasn't able to configure it to work with my own PostGIS Tables. The processing and queries assume that the data is from osm.

Currently I am looking into PGRestAPI, I have been able to create my tiles and host them pretty easily but I am not sure how I will render them on the mobile. I am not able to parse the vector tile endpoints that I get from the server. The Leaflet Vector-Tile plugin is also something I have looked into but I'm not sure how to go about using it on mobile.

How do I go about rendering the vector tiles on a mobile? And is there any alternative way to go about creating a vector tile base map for use on a mobile?


texts: Updating location via data stream using Leaflet

texts:
I'm developing a application where I need to real-time update a location of a map using a data stream. I'm planning to use web-sockets and Leaflet.
Are there any guides I can follow or can you direct me in the correct path?


texts:
How do I select points by expression and delete these selected points?

Is it possible to do this inside the modeler with a tool?

The selection by expression works fine but do I have to set all the selected points to NULL?


texts: Read latitude/longitude to create new point grid at higher resolution

texts:
New to SE, needed to reformat my other question. Basically I'm looking at reading a CSV with longitude and latitude points (in this case housing exposure, simplified below):

HUC4,Longitude,Latitude
1204,-95.956,30.05
1204,-95.955,30.048

And I would like to create a grid for each point in the form below. In my case I need to take the cell resolution (roughly 0.001) and increase this to around 30m or roughly 0.00025 decimal degrees.

newres = 0.00025
csvinput = open(file, 'r')
reader = csv.DictReader(csvinput, delimiter=',')

dst_file = open("out.csv",'wt')
dst_file.write( "longitude,latitude,huc4\n" )
lformat = '%.5f,%.5f,%5s'

for row in reader:
    huc4 = row['HUC4']
    # Longitude Value
    x = float(row["Longitude"])
    xup = x + newres
    xdwn = x - newres
    # Latitude values
    y = float(row["Latitude"])
    yup = y + newres
    ydwn = y - newres

xup  ,ydwn | xup, y  | xup , yup
x    ,ydwn |  x,  y  | x   , yup
xdwn ,ydwn | xdwn,y  | xdwn, yup

Then I just append each row to a new CSV. It's not very efficient.

I'm curious if there's an easier way to do this that doesn't involve raster creation and re-sampling. Ideally I would think that I can take the x/y values and create an array that I can just multiply by

longarray = [xup, x, xdwn, xup, x, xdwn, xup, x, xdwn]
latarray = [ydwn, ydwn, ydwn, y, y, y, yup, yup, yup]

Bottom line, I need to make sure that I am increasing the resolution of these points by a certain factor and I'm scared that if I create a grid out of the points and use GDAL, I will get some floating point errors later.

Any help / pointing me towards the right place?


texts: Why is a GeoTIFF uploaded to MapBox as a tileset shown in black and white and not the original colors in the TIFF?

texts:
When uploading a GeoTIFF file to Mapbox to create a tileset, it is being displayed in Black and White? I would expect it to have the same colors as the TIFF when viewed in any other viewer?

How do you create a tileset with a GeoTIFF to use it as a base map in MapBox?

My ultimate goal is to use this GeoTIFF and then dynamically lay lat/long points on it, but I need to get the map looking correct first.

Seems to work just fine in QGIS.


texts: How to plot pies or histograms around points (GeoPandas) containing different level values indicators?

texts:
This is my code (a simple sample of df). I'd like to have around each point on a map, an histogram or pie showing the mechanical bikes out of the total.

import pandas as pd

pd.DataFrame(
    {'station_code' : ['13014', '13011' , '13012'],
     'numbikesavailable' : ['25', '15', '20'],
     'mechanical' : ['10', '5' , '8'],
     'geometry' :
         ['POINT (2.36404 48.84252)',
          'POINT (2.35822 48.83541)' ,
          'POINT (2.40932 48.8691)']
    }
  )

How can I do this?


texts:
I want to use R in QGIS but I'm unable to find any sign of "R" in Processing  > Processing Toolbox > Options.
How can I install or maybe add R in QGIS? I dont know what this process is called, I'm new in R and QGIS software.


texts: How do I publish a custom .NET GP Tool as a GP Service on ArcGIS Server 10.1?

texts:
I have a GP Tool that I am able to publish as a GP service at 10.0.

I recompiled the DLL for 10.1 and have followed the steps in the Deploying custom .NET and C++ tools page of the 10.1 help doc.  The tool works correctly in arcmap.

However, when I get to the step where I create the geoprocessing package, or when I try to publish directly as a service, I get error 00151 "Tool is unsupported on server".

Has anyone been able to take a custom GP Tool developed with C# and publish it to ArcGIS server 10.1 as a GP service without wrapping it with python?

(I do see this thread about difficulties with an executable - which I presume means an exe.  I would like to avoid having to wrap all my tools in python, as suggested in that thread.)

Update

I probably should mention I also have desktop installed on the same machine with server (along with Visual Studio).  It is interesting that even though I get the 00151 error, when I try to recompile the dll, Visual Studio fails, saying it can't delete the dll because it's locked by another process.  It turns out that other process is arcgis server.  If I stop the arcgis server instance I can recompile, implying that even though I get the 00151 error, arcgis server is still loading the dll into memory.

Update 2

Esri was able to provide me a simplified version of my gp tool.  See the Deploying a Custom C# GP-Tool on 10.1 Server thread in the Esri forum.  The project team will be deciding whether or not to port to python before I spend more time on trying to figure this out.  Will update again once I get access to the test server.


texts:
I have a weird error when using r.watershed in QGIS and I have depleted all my ideas.

Yesterday it was working normally and generating all the files, then suddenly it stopped working and shows the message:

  The following layers were not correctly created:
   Number of cells which drain to each cell
   Direction of the drainage and so on...

I have tried reinstalling the application, re-installing a new version, using GRASS 7 instead of 6, using other DEM, changing the location of the DEM, changing the reference system of the DEM, and the result remains the same.

In addition I should mention that other GRASS algorithms work normally and I get outputs (I had not tried all of them).

Also, before the problem presented, the application started 'mini-dumping' when closing.


texts:
I have the Trial Version for ArcEditor 10.0 with the extensions. Using Maplex/Labeling Manager, after clicking Property, I go to "FITTING STRATEGY". all the lables (stack, overrun, reduce font size, etc..) are all ghosted. I am not given an option to check/uncheck or make changes to anything on that page, nor LABEL POSITION / CONFLICT RESOLUTION pages.

While in Trial, are you restricted what you can do? I thought you were able to use the full Maplex extension, or am I mistaken?

What are my options, if any to make these areas workable?


texts: Closing shapefile using org.geotools.data.shapefile.ShapefileFeatureReader gives warning?

texts:
When working with shapefiles, i get this "warning"

    WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.ShapefileFeatureReader finalize
WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.ShapefileFeatureReader finalize
WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.ShapefileFeatureReader finalize
WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.ShapefileFeatureReader finalize
WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.ShapefileFeatureReader finalize
WARNING: There is code leaving shapefile readers unclosed, this might result in file system locks not being cleared. File is: chranene_uzemi_cr
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.files.ShpFiles logCurrentLockers
SEVERE: The following locker still has a lock: read on file:/F:/Škola/test_data/lesy_cr.shp by org.geotools.data.shapefile.shp.ShapefileReader
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.files.ShpFiles logCurrentLockers
SEVERE: The following locker still has a lock: read on file:/F:/Škola/test_data/lesy_cr.shx by org.geotools.data.shapefile.shp.IndexFile
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.files.ShpFiles logCurrentLockers
SEVERE: The following locker still has a lock: read on file:/F:/Škola/test_data/lesy_cr.dbf by org.geotools.data.shapefile.dbf.DbaseFileReader
Dub 08, 2018 11:23:13 ODP. org.geotools.data.shapefile.files.ShpFiles logCurrentLockers

But the project is build after that and i get results what i need. Not sure where is the problem. I am disposing datastore at the end but it dosnt seems to work.

String overlayPolygons(String shapeFileURL) throws IOException {

String areas = "Object : Area of overlay";

ShapefileDataStore sfds;
sfds = new ShapefileDataStore(new URL("file:///F:\\Škola\\test_data\\lesy_cr.shp"));

SimpleFeatureSource fs;
fs = sfds.getFeatureSource("lesy_cr");

ShapefileDataStore sfds2;
sfds2 = new ShapefileDataStore(new URL("file:///F:\\Škola\\test_data\\chranene_uzemi_cr.shp"));

SimpleFeatureSource fs2;
fs2 = sfds2.getFeatureSource("chranene_uzemi_cr");

SimpleFeatureIterator sfi = fs.getFeatures().features();
double sum = 0;
while (sfi.hasNext()) {
    SimpleFeature sf = sfi.next();
    MultiPolygon mp2 = (MultiPolygon) sf.getDefaultGeometry();
    Polygon p2 = (Polygon) mp2.getGeometryN(0);

    SimpleFeatureIterator sfi2 = fs2.getFeatures().features();
    while (sfi2.hasNext()) {
        SimpleFeature sf2 = sfi2.next();
        MultiPolygon mp3 = (MultiPolygon) sf2.getDefaultGeometry();
        Polygon p3 = (Polygon) mp3.getGeometryN(0);
        Geometry p4 = p2.intersection(p3);
        if (p4.getArea() != 0) {
            sum += p4.getArea();
            areas = areas + "\n" + p4.getArea()+ " : " + p2.getArea() + " : " + p3.getArea();
        }

    }
}

sfds.dispose();
sfds2.dispose();

return "Objects found: " + areas + "\nTotal sum: " + sum;
}


texts:
I am new to OL and I'm trying what I thought was a pretty simple thing. I would like to see a second layer that shows lat/lon points.  Is the code below the correct way to achieve it, or am I missing something obvious?

const projection = getProjection('EPSG:4326');
const view = new View({ projection });
const mySource = new VectorSource({ wrapX: false });
const myLayer = new VectorLayer({ source: mySource, visible: true });
this.map = new Map({target: 'myMap', layers: [
   this.getMyTileLayer(), // Shows up fine
   myLayer // AFAIK I don't see anything plotted
], view});

// Called externally
addToMap(lat,lon) {
  // I can console.log here and I have verified the lat/lons coming in are correct
  const feature = new Feature(new Point(fromLonLat([lon,lat])));
  this.mySource.addFeature(feature);
}

Update:
  If I add the following code to my addToMap function, I now see my marker, but the original TileLayer has disappeared. I'm assuming these issues are somehow related?

this.map.getView().animate({ center: fromLonLat([lon, lat]), zoom: 10});


texts: How to pass arguments to query function in different module

texts:
I have a query (see this jsbin example) to find the percentile rank of a fields value. So it finds all records (N) and then how many records are less than itself (k). Percentile = k/N as a rounded percent.

The math is working correctly, but I would like to use the whole thing a function to which I could pass arguments for field name and field value when you click on features.

My problem is this: I have the listeners for selecting features outside of the AMD module where I have the query. So if the listeners call this percentile function, it is undefined. However, if I move this function outside the AMD module, then query and QueryTask are undefined, because I set those up in the module.

I'm new to Javascript modules, so I'm hoping to find out if there is a simple way I could avoid this problem. I've tried moving the listeners inside the AMD module, but that makes other functions undefined.

I may not be using the word module correctly here. I think what I'm experiencing is a scope problem. I have declared this function in one scope and I'm trying to access it in another. Is that the right way to think about it?

Here is the complete code:

function pctile(pctileField,pctileValue) {

    queryTaskN = new QueryTask("http://gis.datadrivendetroit.org/arcgis/rest/services/DetroitNeighborhoodMetrics/FeatureServer/0");
    queryN = new Query();
    queryN.where = '1=1';
    queryTaskN.executeForCount(queryN, storeResultsN);

    function storeResultsN(results) {
        N = results;
        console.log("N:",N);
    };

    queryTaskN.on("execute-for-count-complete", function() {
        queryTaskK = new QueryTask("http://gis.datadrivendetroit.org/arcgis/rest/services/DetroitNeighborhoodMetrics/FeatureServer/0");
        queryK = new Query();
        queryK.where = '"'+pctileField+'" > ' + pctileValue;
        queryTaskK.executeForCount(queryK, storeResultsK);

        function storeResultsK(results) {
            k = results;
            console.log("k:",k);
            p = k / N;
            console.log(p);
            P = p*100;
        };
    });
};


texts: Sum of point attributes inside hexagons using QGIS UI

texts:
I have a point shapefile that corresponds to locations of animals captured. I've joined a table to this shapefile with additional information on each attribute, namely nb. animals caught in each point. I have another shapefile (hexagonal grid) and my goal is to sum the nb. animals caught for each point that falls inside each hexagon. In a preliminary search on the web, I found that maybe "join by location" could help me do this, but the problem is that I'm getting an "incorrect field name" error, for a field that doesn't even exist in any of the original tables and, on the other hand, it won't let me select the specific field I want it to sum..

Any suggestions? Is there any way I can do this without going into writing scripts mode?


texts: Unable to activate analytic asset in satellite Imagery

texts:
I was trying to import the map of Melbourne. I was able to fetch the details from API coordinates but when i am trying to check the analytic asset status, it is returning inactive.

# extract image IDs only
image_ids = [feature['id'] for feature in search_result.json() .
['features']]
print(image_ids)

# For demo purposes, just grab the first image ID
id0 = image_ids[0]
id0_url = 'https://api.planet.com/data/v1/item-
types/{}/items/{}/assets'.format(item_type, id0)

# Returns JSON metadata for assets in this ID. Learn more:
planet.com/docs/reference/data-api/items-assets/#asset
result = \
requests.get(
id0_url,
auth=('APIKEY', '')
)

# List of asset types available for this particular satellite image
print(result.json().keys())

# This is "inactive" if the "analytic" asset has not yet been
activated; otherwise 'active'
print(result.json()['analytic_dn']['status'])

And That gives me:

['20180831_235330_1008', '20180831_000115_0f31']
dict_keys(['analytic', 'analytic_dn', 'analytic_dn_xml', 'analytic_sr',
'analytic_xml', 'basic_analytic', 'basic_analytic_dn',
'basic_analytic_dn_nitf', 'basic_analytic_dn_rpc',
'basic_analytic_dn_rpc_nitf', 'basic_analytic_dn_xml',
'basic_analytic_dn_xml_nitf', 'basic_analytic_nitf',
'basic_analytic_rpc', 'basic_analytic_rpc_nitf', 'basic_analytic_xml',
'basic_analytic_xml_nitf', 'basic_udm', 'basic_udm2', 'udm', 'udm2'])
inactive

The status is returning as inactive. It should be active.


texts: Which character encoding is used by the DBF file in shapefiles?

texts: Which character-encoding is used by the dbf-file in shapefiles? It seems it is handled different, based on the program and the local encoding-settings of the machine. Which encoding is 'right' - specified for the format?


texts: ST_Distance_Spheroid returning very large value

texts:
I have hard coded here one of the coordinates. In reality, the distance in km between these two points is around 2.2 km but the query returns 8033830.232829098 which as suggested by the PostGIS documentation should represent the distance in meters => 8033 km is a lot larger.

select ST_DistanceSpheroid(
            ST_GeomFromGeoJSON(
                '{"type":"Point","coordinates":[53.4623013676122,-2.25182451244364]}'
            ),
            ST_GeomFromGeoJSON(
                (
                    select
                        coordinates::text
                    from
                        actor
                    where
                        actor.id = '4b9722c0-50dd-11e7-b8e3-4d48263ab9fd'
                )
            ),
            'SPHEROID["WGS 84",6378137,298.257223563]'
        )

What do I miss ? I have checked various threads online and the query seems fine to me.


texts: How to create ArcGIS Online map to only only show intersecting line features when clicking on a point feature

texts: I have a map with a feature class containing airport locations and another feature class showing flight paths between airports. I would like to create a map on ArcGIS Online where when I click on a particular airport, only flight paths from that specific airport are shown. Is there a way of doing this in ArcGIS Online ideally without too much coding? This is a sample of the map in it's present state. I would like to only show the flight routes from specific airports when they are clicked, for example, if I click London, only show flight paths from London. I was having a look through the ArcGIS Online sample code but couldn't find anything with this kind of functionality.


texts:
I have 2 raster maps both with the same size and projection etc.
One map is some kind of base map (jpeg format) and the second map is an overlay that adds something (border lines, points...) to the base map (png format with alpha channel). I want to merge both maps together, i.e. the second map shall overlay the base map with a transparency dependent from its the alpha channel. ImageMagick composite does exactly this but the geo information get lost. Is there any option to do the job with gdal tools?
I tried already gdalmerge and creating a vrt file, but in both cases just one map was copied to the output.


texts: GDAL, Python, pyramids/statistics read/write

texts:
As I work through learning programing with raster imagery, I am trying to read in a raster, change the pixel values and write the output file. Below is the code. I have blocked out the change pixel values as my output file was empty and in the wrong location, so I am now trying to just read in a file by block and write out the same file. However  the file I get is defaulted to the 64/64 block size, and my skip factor is 999, so my output raster does not 'match' my input raster. I have not been able to SET the metadata changes. Can someone tell me the correct way in a script to set new metadata options, and then I will need the numpy.where to process correctly?

import os, sys, ogr, gdal, utils, numpy
from gdalconst import *

file = (r'C:\Users\Desktop\images\large_raster.img')

gdal.AllRegister()

ds = gdal.Open(file, GA_ReadOnly)
if ds is None:
    print 'Could not open file'
    sys.exit(1)

rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

band = ds.GetRasterBand(1)

yBlocksize = 512
xBlocksize = 512
count = 0

driver = gdal.GetDriverByName('HFA')
metadata = driver.GetMetadata()

#print metadata

outfile = (r'C:\Users\aschilli\Desktop\images\chrisG_test2.img')

outDataset = driver.Create(outfile, cols, rows, 1)

outBand = outDataset.GetRasterBand(1)
outBand.FlushCache()
outBand.SetNoDataValue(-99) #this gives 0
outDataset.SetGeoTransform(ds.GetGeoTransform())
outDataset.SetProjection(ds.GetProjection())
gdalinfo = (r"C:\Python27\ArcGIS10.3\Scripts\new_ve_folder\Lib\site-packages\osgeo\gdalinfo.exe")
command = gdalinfo + ' -hist ' + ' -stats ' + outfile
os.system(command)
gdaladdo = (r"C:\Python27\ArcGIS10.3\Scripts\new_ve_folder\Lib\site-packages\osgeo\gdaladdo.exe")
command2 = gdaladdo + ' --config ' ' HFA_USE_RRD ' ' YES ' + outfile + ' 3 9 27 81'
os.system(command2)
#this code below does nothing
outDataset.SetMetadataItem({
    'BLOCKSIZE': '512/512',
    'STATISTICS': 'TRUE',
    'DMD_CREATIONDATATYPES': 'Int16'
})

#loop through rows
for i in range(0, rows, yBlocksize):
    if i + yBlocksize < rows:
        numRows = yBlocksize
    else:
        numRows = rows - i
#loop through coloumns
    for j in range(0, cols, xBlocksize):
        if j + xBlocksize < cols:
            numCols = xBlocksize
        else:
            numCols = cols - j

    #read data and do calculations
        data = band.ReadAsArray(j, i,  numCols, numRows)
        # out_data = numpy.where(abs(data) < 500, 10, data)

        outBand.WriteArray(data, j, i)

del outDataset

#http://gis.stackexchange.com/questions/160578/write-python-ndarray-to-raster?'


texts: Creating cross-sections in ArcGIS Desktop?

texts:
Is there the possibility if I have a polyline (named "Linedivide") which is divided into,  let's say, 5 parts to make ArcGIS draw perpendicular lines of special length at the middle of all 5 parts of the Linedivide?

Can someone tell me where I can find a tool for this in the toolboxes?


texts: Problems with batch processing a raster calculator

texts:
I have a lot of (75) .asc files which I want to do the same simple raster calculation on each of (in QGIS). My calculation works when I use the raster calculator on one layer/file, but I am struggling with how to get it to run as a batch process across all 75.
When I open the raster calculator from the processing toolbox and click 'run as batch process', I assume I need to be entering my formula in the 'expression string builder'. Is that right?
If so, I don't see how to enter simple expressions (e.g. log10 of the value). If I should be doing the calculation on the 'raster_value', (a) what should I enter for its 'point' argument? (b) how do I avoid having to type the layer name manually for each of the 75 files and (c) how do I get it to set a the output name automatically (i.e. not have to manually type a different one 75 times)?
I've tried reading the manual/searching here without success.
Or would it be better to do this outside of QGIS?
Thanks very much!
John
[EDIT - added screenshots to show what I'm talking about. First is batch process window of raster calculator; second is what I assume to be the way to add an expression; third is the window where I assume you enter the expression, but as per the question above, this is where I'm stuck. Thanks!]


texts: Missing dependency of grass7.2, QGIS2.18 ;r.mapcalculator

texts:
I install QGIS2.18 with GRASS7.2.0. QGIS worked well until I was about to used the 'r.mapcalculator' through grass.

I am unable to find it though the 'processing' tab, and type the program states that the plugin is a missing dependency "GRASS folder is not configured. Please configure it before running GRASS algorithms."

How do I configure it?


texts: What is the better way to improve performance of a geocoder, faster disk or more RAM

texts:
I'm using free tiger geocoder for getting coordinates, but the performance is really slow for the bigger states.
E.G.

A query for a small state can take 300ms while a query for a big state can take about 1min.
My guess here is that the machine runs out of memory for big states.


texts: ValueError geopandas.sjoin: why does it not read my GeoDataFrame?

texts:
I'm working on a spatial join between a polygon from a shapefile and points from a CSV file, using geopandas.sjoin. I want to assign each point with a name of the province it is located in. I do know how to do it in R, but I'm restricted to use Python in this context.
I found an example of a sjoin over here, so I tried it:

My data looks as follows:

import pandas as pd
import geopandas as gp
from shapely.geometry import Point
from geopandas import GeoDataFrame

#My observation points:
obs = pd.read_csv('data/butterflies.csv')
pts = gp.GeoDataFrame(obs, geometry=geopandas.points_from_xy(stations.lon, stations.lat))

pts.head()

#Polygons of some provinces in The Netherlands:
nlpoly = gp.read_file('data/provinces_NL.shp')
nlpoly

Now when I try the spatial join as was suggested in the example above, it returns me a ValueError:

gp.sjoin(pts, nlpoly[['name']],how='left',op='intersects')

>>ValueError: 'right_df' should be GeoDataFrame, got <class 'pandas.core.frame.DataFrame'>

I checked the types for both objects but they were all GeoDataFrames as it should be, yet sjoin does not recognize my rigth_df as a GeoDataFrame.
I guess the bug is in the referencing of the column "name" that I want to assign to the observation points, but I have no clue how to do it right.


texts: Count number of pixels in range of values in GEE

texts:
I have raster with NDVI values and I would like to use reducer in order to know how many pixels in the raster are smaller than 0.2, between 0.2-0.4 , betwen 0.4-0.6 and gt than 0.6.

i'm not sure that the way I have been that gives me true results, because I have difference between the total number of pixels to the sum of the differnt categories pixel counts.

    var c0406 = image.gt(0.4).lt(0.6).selfMask().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 20,
  }).values().get(0);

when each time I have change the condition in image.gt

I'm not sure it's correct because when I calculte the total number of pixels, I get 17,748 , but when I calculte the categories pixels I get around 45k:

   var total = image.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 30,
  }).values().get(0);

*i'm not sure what is systen asset size and if it's also suppose to be number of pixels

My end goal is to count correctly the number of pixels between different ranges (0.2-0.4 ect.)


texts: Can privileges be set to allow certain users read/write access to a protected version?

texts: I have multiple versions of a database (some protected and some public). Would it be possible to allow some of the database users read/write access to the protected versions? I have read of protected versions; Any user may view the version, but only the owner or the geodatabase administrator may edit datasets in the version or the version itself.


texts:
I want to make something like in this question in QGIS.

In ArcGIS there is available Make Divided Roads plugin. Is there easy way to make the same job in QGIS?


texts: Looping Summary Statistics on changing field name

texts:
I am trying to calculate summary statistics for a number of shapefiles, which is easy enough, but the field name for each shapefile is different. The field name is identical to the file name if that makes things easier.

I'm trying to use the following script:

import arcpy, os
from arcpy.sa import *

arcpy.env.workspace = r'C:\Users\kellyj\Desktop\Projects\CostAnalysis\SHP'

for file in arcpy.ListFiles("*.shp"):
    arcpy.env.overwriteOutput = 1
    filename = file[0:-4]
    print "Calculating " + filename
    arcpy.Statistics_analysis(file, arcpy.env.workspace + "/" + filename, [[filename, "SUM"]], "")

print "Done"

How can I write this script so that I can use the file name as the statistics field?


texts:
I try to export a GML file from HALE editor but the validation has always errors with this text:

 The value '' of attribute 'gml:id' on element 'ps:ProtectedSite' is not valid with respect to its type, 'ID'., on line_, column_

  '' is not valid value for 'NCName' on line_, column_

Anybody can help me how to fix this errors?


texts:
I got an exception when exporting a .tif image from QGIS.

Cannot write raster. Error code: data source creation
Cannot create new dataset  mymap3.tif:
Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

What can the problem be?


texts: String Formatting;-Calling values from a dictionary to select by attribute in fc

texts:
New to python. In the code below, I create a dictionary b={} and wish to call keys to select by attribute. All goes well until I am required to format strings to come up with whereclause. I am constantly getting an error.

 b = {}
    with arcpy.da.SearchCursor(tc, (field)) as rows:
        for row in rows:
            if row[0] not in b: #I am using 'not in' to test value in dictionary
                b[row[0]] = [row[1],row[2]]

    PATCH_lyr = arcpy.MakeFeatureLayer_management(fc, "PATCH_lyr")
    for k, v in b.iteritems():

            print k, v

   PATCH_lyr = arcpy.MakeFeatureLayer_management(fc, "PATCH_lyr")
for k, v in b.iteritems():

        print k, v

frt = b.iteritems()
for ind in frt:
        whereClause =  "IDO = '{0}'".format(ind)
        kala = arcpy.SelectLayerByAttribute_management(PATCH_lyr, 'NEW_SELECTION', whereClause

)


texts: Does ArcGIS Desktop have equivalent to RStudio Console Windows?

texts:
I know Esri provides scripting capacity, via the Python window, but I've had a hard time with it. Making the jump from the Field calculator to the Python window was too much to me. Trying to keep 3 different sets of syntax straight (Esri, VBA, and Python) was misery (Esri doesn't care about caps or whitespace, Python does...)

Starting use RStudio as a GIS. One of the reasons I'm starting to move from ArcGIS to R is the need for reproducibility in my analysis. Much of the analysis I'm doing is scripted. But wondering if there is a way to keep using Esri.

One of the features I really like RStudio about it is the Console, which acts as a command line interface. But it also logs the equivalent of the command line code you would have typed in when you run something from the GUI. This has been very very helpful in learning R.

Is there an equivalent for the Python window, so I could see the inputs from my GUI in Esri as the equivalent Python output code?

That would make translating GUI inputs into a Python script much much easier. It would also make reproducing analysis much easier, as I'm repeatedly being called on to perform similar analyses, but months apart.

Using 'Record Macro' in EXCEL to record the VBA equivalent of my GUI commands, and then editting the VBA code to create a new script/macro.

Is there a similar way to document my workflow in ArcGIS?

Related:

Developing GUI in Python for ArcGIS geoprocessing using PyQT/Tkinter/wxPython?
Can anyone recommend an Arcpy GUI?


texts: Finding unnamed Antarctic peaks using QGIS?

texts:
I'm a complete novice at GIS. My project is to find the highest peaks in each Antarctic territory, and compare the results with the international gazetteer to find which peaks are and aren't named.

I'm trying QGIS Desktop with Grass with the 200m resolution radar altitude DEM raster file. Locations are (x,y) from the south pole in metres.

ftp://sidads.colorado.edu/pub/DATASETS/nsidc0082_radarsat_dem_v02/200M/ARCINFO/demosu200_v2.tar.gz

The Gazetteer is https://data.aad.gov.au/aadc/gaz/scar/download.cfm and is available as CSV, GeoJSON (whatever that is) and KML formats. Locations are in latitude and longitude to the nearest minute of arc, which is not accurate enough.

Territories are given by longitude range in Wikipedia

I've succeeded in opening the DEM with QGIS. I could find peaks manually one at a time using colours and the "identify features" button.

But is there a better way?


texts: Get all geom points that are more than 3 meters from the linestring, at big scale

texts:
With the help of Osmosis i was able to import all the data of a Open Street Map city, (Nodes, ways, etc...).
In another table at same db, i have Millions locations in that city, these locations are also geom points.
The image below shows a way that is a linestring and the green dots represents locations that are geom points.

Is it possible to query at really big scale without iterating over each location, to get all the locations (geom points) that are more than 3 meters away from its closest way (linestring) within a all the ways i have in the db?
My first approach was this query
SELECT locations.geom FROM locations
JOIN ways
ON ST_DWithin(
        (ways.linestring)::geography,
        (locations.geom)::geography,
3) = false;

and
SELECT locations.geom FROM locations,ways
WHERE ST_DWithin(
    (ways.linestring)::geography,
    (locations.geom)::geography,
    3) = false;

with these two indexes:
CREATE INDEX ways_geography_cast_gpx ON ways USING GIST (geography(ways.linestring));
CREATE INDEX locations_geography_cast_gpx ON locations USING GIST (geography(locations.geom));

But it takes lots of time.


texts: Joining tables to feature classes in batch using ArcPy?

texts:
I need some help with batch processing using python. I am trying to join fields to multiple point features classes from multiple dbf tables (650 feature classes and tables). Names of the dbf tables match names from each feature class. I have got the script so far; I am not sure where I went wrong as the script runs but feature classes are not updated with the join. I am using arcpy.joinField.management instead of Add Join because I want to update the feature classes directly without having to create feature layer and save.

Is that where I am wrong?

My script so far is below:

 import arcpy, os
 from arcpy import env
 # Add variables
 GDB1 = r"C:\test\GDB1.gdb"
 GDB2 = r"C:\test\GDB2.gdb"
 joinField = "Work_ID"
 env.workspace = GDB1
 #Join tables
 for fc in arcpy.ListFeatureClasses():
     env.workspace = GDB2
     for table in arcpy.ListTables():
         try:
             if table() == fc():
                 try:

                        arcpy.JoinField_management(GDB1 + os.sep + fc,      joinField, GDB2 + os.sep + table, joinField, "")

                except:
                       arcpy.AddMessage("Could not join ")
                       pass
    except:
            pass
env.workspace = GDB1


texts: Loading shapefile into CartoDB throws 2001 error?

texts:
I've been trying to load shapefiles (that are zipped in windows compressed folder) into CartoDB and I always get an error 2001.

These shapefiles were exported out of an ESRI geodatabase (SQL Server).  I believe they are in the correct projection - Web Mercator auxilliary sphere.

Has anyone had a similar experience or might know what is happening?


texts: Multiple Labels from attribute column in QGIS

texts:
I was wondering if there is a way to have multiple labels from a joined attribute table in QGIS. I have two tables, a polygon table, and a table with words (codes). The codes table can be joined by a foreign key (geoid) with the codes table. I would like to simulate a word cloud (see image below) in QGIS doing the following (if it is possible):

1) Set a label from a joined table
2) Have many labels shown, and sized according to their frequency in the table

Simulated word cloud

Codes table with words to display and frequency of words.

Any thoughts or insight is welcomed.


texts: Calculating the area of groups based on intersect of several shapefiles using ArcPy?

texts:
I have 5 shapefiles. There are 72 possible groups of 'land cover+top-rock+subrock+soil drainge' in each watershed. I want to get the area of the different groups of 'land cover+top-rock+subrock+soil drainge' in each watershed.

I have thought about:

Intersect the shapefiles.
Create a new column named group ='land cover+top-rock+subrock+soil
drainge' columns
Calculate area of each polygon in the column 'group'

Is there a way to get the area of groups resulting from intersect of several shapefiles at once using ArcPy?


texts: Creating buffer only in specific direction using ArcGIS for Desktop?

texts:
I am trying to create a buffer for several polygons in a south-western orientation. As far as I know, this is not possible using the buffer tool (I use ArcGIS 10.3). I could do it manually but for some 400+ polygons it would take far too long.

Does anybody know a better way?

This is more or less what I am aiming for:


texts: Identifying ArcMap document version without ArcPy?

texts:
Is it possible to identify what version (10.0, 10.5 etc.) an ArcMap document (.mxd) has been saved in?

I am aware of ArcPy method to determine ArcMap document version, but I am interested in a solution that does not require ArcPy.

I was hoping for a solution that was easy enough for any user to perform; a "right-click-and-properties" kind of solution. I thought something like this would be a given feature in ArcCatalog, but I cannot seem to find the information there nor in Windows explorer.


texts: Assign a point to polygon using pandas and shapely

texts:
I have a series of Easting and Northing coordinates in pandas dataframe along with other columns of information. I know the coordinates fall within and administrative area and want to match the coordinates and rest of columns to that administrative area, I have the shapefile of the administrative boundaries. I'm familiar with R and visually this looks easy, I can plot the coordinate and see it fall into the area. However i have n=8,335 records to match and need to use Python as all other records are loaded and the analysis i have is better suited to Python.

So far i have work flow as;

Convert Easting and Northing to lat and long

Somehow use Shapely to tell me if point is in polygon / I understand that matplotlib.path.Path.contains_points executes faster.

Assign point to polygon and add new column to dataframe with the polgon reference in it.

Use pandas merge to join dataframe to others. Project to map again and check visually.

Does that sound like a sensible approach or is there a more efficient route?


texts: QGIS Georeferencer plugin TIF format problems

texts:
I'm recently new to QGIS, and have been using it to georeference some photos for a physics thing. Sadly I seem to have a pretty poor handle on using TIF files with it.

When I load up a jpg image, it looks the same as when I load it up in say paint, but when I try to georef a TIF image it becomes very discolored and more difficult to use.

I would like to know if there is any way around this such that I can use the original TIF image to georeference. Right now, I convert it into jpg format, georeference it, save my points and then load up the weird looking TIF format with those points to get a georeferenced TIF image. I'm worried that I'm losing data when doing this though due to the weird colors of the image. It's very important that I keep as much data as possible, hence why I can't just convert to jpg anyways.


texts: Bulk grant privileges to multiple tables at once in ArcGIS Desktop?

texts:
I am able to grant multiple privileges (SELECT, INSERT, UPDATE, DELETE) to multiple Users/Roles to one single Table each time on an Enterprise Geodatabase from ArcCatalog 10.1.

That's done by Manage -> Privileges... for each Table.

Now that I have got a large number of Tables/Views, just wondering whether there's a way/tool/PyScript of bulk granting privileges to all the Tables/Views at once.


texts: How to implement a simple QFileDialog in a QGIS plugin

texts:
I'm building a QGIS plugin which only requires user interaction to select a save location for an output file. I used Plugin Builder to create a template and assumed I could simply replace the template dialog with a QFileDialog object. Here's the relevant code:

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

    # Create the dialog with elements (after translation) and keep reference
    # Only create GUI ONCE in callback, so that it will only load when the plugin is started
    if self.first_start == True:
        self.first_start = False
        self.dlg = QFileDialog(self.iface)
        self.dlg.setFileMode(QFileDialog.AnyFile)
        self.dlg.setNameFilter(self.tr("Zip files (*.zip)"))
        # self.dlg = GiftWrapDialog()

    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        pass

This code throws an error when I try to run the plugin:

TypeError: arguments did not match any overloaded call:   QFileDialog(QWidget, Union[Qt.WindowFlags, Qt.WindowType]): argument 1 has unexpected type 'QgisInterface'   QFileDialog(parent: QWidget = None, caption: str = '', directory: str = '', filter: str = ''): argument 1 has unexpected type 'QgisInterface'

Is there any way for my plugin to simply show a standard file dialog without having to create my own dialog?


texts: Generate clusters of X points within N meters

texts:
I have a set of roughly 4,000 points in very close proximity and I need to generate clusters of 5 points, with each point no farther than 50 meters from the (mean?) center of the cluster. I've explored various tools in ArcMap, especially in the analysis and spatial statistics toolboxes, though I can't seem to find a tool that will do the trick. I know that I would need roughly 800 clusters, but I'm not sure how to derive them from the point data.

The buffer solution mentioned here won't work because the points in my data set are in such close proximity that buffering each point and then dissolving would create one or two large polygons that cover most of the data set. Also, I need to each point in my data set to be assigned to a cluster, while creating as few clusters as possible (a small number of clusters with less than five points is acceptable).


texts: GPL-compatible database of country & region names

texts:
For WordPress's sister project BuddyPress, we are looking for relatively-reputable resources that manage a list of countries and regions. Due to the WordPress OpenSource philosophy, we are only seeking GPL-compatible resources.

If no such resources exist, we could possibly consider GPL-compatible services as well.

Or there could possibly be a GlotPress-like service (see http://translate.wordpress.org/projects/wp/dev) wherein the community contributes and maintains the database which can be downloaded via a link or API.

So far all the databases/services I have found have either been commercial or otherwise incompatible with GPL.


texts: Combining field data from multiple shapefiles into one queryable layer using ArcGIS for Desktop?

texts: How do you combine multiple field collected data set in shapefile form in to a single layer in ArcGIS 10.3 for Desktop?


texts: What kind of file formats supports the torque marker-file?

texts:
I want my marker-file to be a gif 'cause I want to have a moving footsteps effect.

I've tried but it doesn't work well, it shows like the gif is an static image.

  tami_3[value=1]{marker-file: url(https://thumbs.gfycat.com/GratefulPlainBactrian-size_restricted.gif);
  }

When I try with other gifs, they move, but they are really small.
Someone with the same problem?


texts: Map two sf objects on one map with mapview

texts:
I have two sf objects, precincts and elemZones. I want to make a mapview object showing the boundaries of both sf objects, and I want to style each sf object differently. My code so far.
library(sf)
library(mapview)

precincts <- read_sf("precincts.shp")
elemZones <- read_sf("elemZones.shp")

mapview(
  x = elemZones,
  legend = F,
  color = "black",
  alpha.regions = 0,
  popup = F
)

This only maps elemZones. How do I add precincts?


texts: Importing massive amounts of raster files

texts:
I have 17.000 individual files that I would like to load into QGIS.

These files are .img-files (Rasters) that contain land elevation data. Each file is about 1 Mb in size. I tried to add in all the layers, but QGIS would not respond.

How can this data be imported to QGIS?

I'm using QGIS 3.4.10


texts:
I'm developing a C++ application and I have to parse a OGC Filter in order to implement querying functionality on a metadata repository.
Does anyone know any C++ library to parse OGC Filters?

Thanks.


texts:
I have this shell script:

#!/bin/sh
r.mapcalc water = 120
d.rast water
i=1
while [ $i != 100 ]
do
    n=1
    while [ $n != 10 ]
        do
            r.mapcalc < water.mapcalc
            d.rast water
            n='expr $n + 1'
        done
        g.copy rast=water,water.$i
        i = 'expr $i + 1'
    done

It is from http://grass.osgeo.org/gdp/raster/mapcalc-algebra.pdf. I understand the syntax of the code since I am a programmer, but I cannot understand the line r.mapcalc < water.mapcalc. What does it mean? Sorry first time to use a shell script.


texts: How to find the second highest value in a rasterStack in R

texts:
I have a rasterStack of rasters representing count data. I want to return a single raster whose values are the index of the raster with the second highest value for each cell. (I used which.max for the highest value.)

I'm using the raster package in R, but could use ArcGIS too. Thanks!

This works for identifying the raster with the max value:

set.seed(123)
library(raster)
r1 <- raster(nrows=10, ncols=10)
r2=r3=r4=r1
r1[]= runif(ncell(r1))
r2[]= runif(ncell(r1))+0.2
r3[]= runif(ncell(r1))-0.2
r4[]= runif(ncell(r1))+0.5

rs=stack(r1,r2,r3,r4)
plot(rs)

rs.max<-which.max(rs)
plot(rs.max)


texts:
Is it possible create a cluster source with tiles?

With geojson works perfectly but I would like to work with a tile source.

map.on('load', function () {

   // add source
   map.addSource('geosjonsource', {
    'type': 'geojson',
    'data': 'https://docs.mapbox.com/mapbox-gl-js/assets/earthquakes.geojson',
    'cluster': true,
    'clusterMaxZoom': 14,
    'clusterRadius': 50,
   });

   // add layer
   map.addLayer({
     'id': 'prueba',
     'source': 'geosjonsource',
     'type': 'circle',
     'filter': ['has', 'point_count'],
     'paint': {
       'circle-color': ['step', ['get', 'point_count'], '#51bbd6', 100, '#f1f075', 750, '#f28cb1' ],
       'circle-radius': ['step', ['get', 'point_count'], 20, 100, 30, 750, 40 ]
      },
      'layout': {}
   });
}

I've created a vector tile server to render pbf sources and I would like work with this source:

map.addSource('vectorsource', {
    'type': 'vector',
    'tiles': ['http://192.168.240.71:8090/tile/earthquakes/{z}/{x}/{y}.pbf'],
    'cluster': true,
    'clusterMaxZoom': 14,
    'clusterRadius': 50,
});


texts: Is there a sequence stored in a QGIS layer that can be accessed via pyqgis?

texts:
I have a layer from spatialite. When I add a feature I want to know it's identifier (based on the primary key (integer) which is auto-incremental : for example XXXXXX-72) before the creation on the entity (attributeForm OK + commit) by using pyqgis.

My issue is when i delete one or several features. I can predict the next identifier only if I know how many features I've deleted which will quasi never be the case.

I can only access the number of features or every unique id but that doesn't help much when, for example, the last three features were deleted.

Is it possible to access the sequence of the table that store every unique identifier ?


texts: QGIS Desktop 2.14.0 with GRASS 7.03: Unable to edit

texts:
I have imported two text layers into QGIS but now I am not able to edit them. All the editing fields are grey.
I don’t use QGIS that much, so it’s a stupid question.. But how can I activate the edit functions?


texts: Housenumber(point) angle from line (street)

texts:
using arcgis 10 with Point data (Housenumber without angle attribute) and line data (street without angle attribute)

I want to create housenumber labels from point with correct angle from nearest streetline.

Is there any toolbox for this?.

my thinking of this is calculate AVERAGE streetline angle and then somehow put this values to nearest point?


texts: How to go around QGIS's same CRS restriction when doing spatial queries?

texts:
As someone who used ArcGIS pretty frequently and then moved to QGIS, I find this restriction to be pretty annoying.

Normally it would be possible to convert everything to the same coordinate system but the problem is that I'm working with WFS layers and QGIS goes absolutely bonkers when I try to change the CRS (only way is to ctrl alt del and end task)

I honestly thought that the select by location tool was broken since every selection I've attempted simply would not work.

Is there a work around this?


texts: Connect to QGIS 2.6 Python Plugin function via button on custom feature form

texts:
I am trying to connect to a QGIS 2.6 Python plugin function via a custom feature form button.

On my feature form, I'm attempting to do a few things :

  enable the updateButton only if the checkbox (Reparee) is checked
  disconnect the default validate and reconnect to my validate function
  connect to a function in a python plugin with the click of a button

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis import utils

OriginInfofield = None
Repareefield = None
myDialog = None

myFeature = None
myLayer = None

def formOpen(dialog,layer,feature,featureid):
    global myDialog, outilInstance
    myDialog = dialog
    global Repareefield, OriginInfofield

    myFeature = feature
    myLayer = layer

    OriginInfofield = dialog.findChild(QLineEdit,"OriginInfo")
    Repareefield = dialog.findChild(QCheckBox, "Reparee")
    buttonBox = dialog.findChild(QDialogButtonBox,"buttonBox")
    updateButton = dialog.findChild(QDialogButtonBox, "updateButton")

    # Disconnect the signal that QGIS has wired up for the dialog to the button box.
    buttonBox.accepted.disconnect(myDialog.accept)

    # Wire up our own signals.
    buttonBox.accepted.connect(validate)
    buttonBox.rejected.connect(myDialog.reject)

    updateButton.clicked.connect(update)

    if Repareefield.isChecked():
        updateButton.setEnabled(False)

    if 'fuitesOutil' in utils.plugins:
        outilInstance = utils.plugins['fuitesOutil']

    QMessageBox.information(None, 'Info', 'Got to open form')

def update(featureid):
    global updateButton, Repareefield

    outilInstance.updateDeclaredTableWidget(featureid)

    QMessageBox.information(None, 'Info', 'Got to function')

def validate():
    global outilInstance
    global OriginInfofield

    if not len(OriginInfofield.text()) > 0 or OriginInfofield.text() == 'NULL':
        OriginInfofield.setStyleSheet("background-color: rgba(255, 107, 107, 150);")
        QMessageBox.warning(None,'Probleme de saisie',"Veuillez remplir 'Origin d'Information'")
    else:
        # Return the form as accepted to QGIS.
        myDialog.accept()

My plugin function looks like this (part of it -- the function works fine when triggered within the plugin itself) :

def updateDeclaredTableWidget(self, fid):
        QMessageBox.information(None, 'Info', 'Got to the plugin')
        self.decdlg.show()
        request = QgsFeatureRequest(fid)
        ...

The only QMessageBox I reach is the one in the validate function, but the condition does not work and the dialog is accepted anyways. I'm surprised the one in formOpen doesn't pop up. My python file is in OsGeo4W/apps/qgis/python. Previously, i kept it in .qgis2 but the structure seems different now that i've installed using OsGeo4.

UPDATE BELOW:

I'm looking into the updated QgsAttributeForm class reference and notice that accept and reject are now deprecated. Unless I'm overlooking some other error, I feel like I must be poorly disconnecting the accept button on my feature form. I've tried to replace them with resetValues and Save, but still have had no success.

  Using myDialog.disconnectButtonBox() gives me an AttributeError:    QgsAttributeForm object has no attribute disconnectButtonBox
  I then tried disconnecting the way I had been doing, but with the new save function. This gives me a TypeError: disconnect() failed between 'accepted' and 'save'.

Any pointers on properly disconnecting the buttonBox and reconnecting to my functions?

My attempts so far:

def formOpen(dialog, layer, feature, featureid):
    myDialog = dialog

    #featureForm.disconnectButtonBox()
    buttonBox.accepted.disconnect(myDialog.save)
    buttonBox.accepted.connect(validate)
    buttonBox.rejected.connect(myDialog.resetValues)

def validate():
    if not myConditions == True:
        doStuff
    else:
        myDialog.save()

Still breezing through the condition.

My ultimate goal is to pre-fill some fields based on another intersecting layer's feature attributes.


texts: Change Text Labels to Fit New Size Paper

texts:
I am using ArcGIS 10.2.2 for Desktop.

I am trying to create a new map based on a previous map which was in paper size Arch D. The new map is sized 11 by 17. I have fit the layers into the new maps, but I am unable to change my labels to fit on the new map.

My labels are annotations and right now are at the same scale as they were in the old map.


texts:
I'm attempting to display some properties from the points on my map. I'm using the example provided on the geoserver site: http://openlayers.org/en/latest/examples/getfeatureinfo-image.html?mode=raw. When I click on my map, it opens up a dialog box, but its empty. I checked the logs and nothing seems unusual other than the fact that no matter where I click on the map, pixelX and pixelY are always 50, which seems to be the center point, as the width and height in the request are both 101.

Code

    <title>WMS GetFeatureInfo (Image Layer)</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.2.0/css/ol.css" type="text/css">

&nbsp;

Log

Request: getFeatureInfo
    BaseUrl = http://localhost:80/geoserver/
    Exceptions = application/vnd.ogc.se_xml
    FeatureCount = 1
    Get = false
    GetMapRequest =
GetMap Request
 version: 1.3.0
 output format: image/png
 width height: 101,101
 bbox: SRSEnvelope[-267350.5296874754 : 7638072.6836785935, 2.086165940819164E7 : 2.876708262155771E7]
 layers: cite:points
 styles: customPoint
     InfoFormat = text/html
     PropertyNames = null
     QueryLayers = [org.geoserver.wms.MapLayerInfo@f58b334c]
     RawKvp = {CRS=EPSG:4326, FORMAT=image/png, INFO_FORMAT=text/html, I=50, STYLES=, J=50, WIDTH=101, HEIGHT=101, LAYERS=cite:points, REQUEST=GetFeatureInfo, BBOX=-267350.5296874754,20861659.40819164,7638072.6836785935,28767082.62155771, VERSION=1.3.0, SERVICE=WMS, QUERY_LAYERS=cite:points, TRANSPARENT=true}
     Request = GetFeatureInfo
     RequestCharset = null
     Version = 1.3.0
     XPixel = 50
     YPixel = 50

The fact the the example doesn't work makes me think it has something to do with how I set up the layer or server, but I haven't found anything yet.


texts: Can I use grass r.mask and t.rast.univar in concurrent processing (threads of python?)

texts:
I have this python script:
    if gscript.find_file("MASK", element="cell", mapset=MAPSET)["file"]:
      gscript.run_command("g.remove", flags="if", quiet=True, type="raster", name="MASK")
    gscript.run_command("r.mask", vector=vect, quiet=True, where="gid = "+cat)
    gscript.run_command("t.rast.univar", overwrite=True, input=RASTER_TIME_SERIE+"@"+MAPSET, separator=",", output=output, quiet=True)

But, this processing is long running.
I'd like to execute this routine (for multiple value of "cat") in threads of python. But, I don't know what is the method to run in concurrent threads, multiple MASK values.
Any idea?

PD: I'm not a expert of GRASS. I know python and gis, only

jorge infante
rosario - santa fe - argentina


texts: What is the path for font in mapfish print module for windows in config.yaml?

texts:
The following path is for linux:

fonts:
  - /usr/share/fonts/truetype/ttf-dejavu/DejaVuSans

What is the path in order to call font type from windows?


texts:
I have features in drop down list.when I select feature in drop down it is zooming to that layer.Now I wanted to highlight that zoomed feature.how can I style the zoomed feature?

can any one help me to achieve it?

here is my zooming function:

    var owsURL = "http://192.168.1.137:8080/geoserver/school/ows";
var SearchResultsLayer = new ol.layer.Vector({
    id: 'SearchResultsLayer',
    style: function (feature, resolution) {
        style.getText().setText(resolution < 5000 ? feature.get('') : '');
        return styles;
    }
});

function Zoomto(filter, SearchResultsLayer) {
    var schoolname = $('#ddlschool').val();
    alert(schoolname);
    var filter = "talukname='" + schoolname + "'";
    alert(filter);
    $.ajax({
        url: owsURL,
        type: 'GET',
        contentType: 'application/json',
        dataType: 'jsonp',
        jsonpCallback: "parseResponse",
        data: {
            typeName: "school:taluk",
            outputFormat: "text/javascript",
            service: "WFS",
            version: "1.1.0",
            request: "GetFeature",
            CQL_FILTER: filter

        },
        success: function (response) {

            var SearchResultsLayer = new ol.layer.Vector({});

            SearchResultsLayer.setSource(null);

            var searchresultssource = new ol.source.Vector({
                features: (new ol.format.GeoJSON()).readFeatures(response, {
                    featureProjection: "EPSG:3857"
                })
            });

            SearchResultsLayer.setSource(searchresultssource);
            map.getView().fit(SearchResultsLayer.getSource().getExtent(), map.getSize());

        }
    });
}


texts: Does GPS correct for continental drift, and if so, how? (I'm mainly interested from the perspective of consumer devices such as smartphones.)


texts: Queries about Optical Remote sensing images

texts:
I am new in Remote sensing and studying about Optical Remote sensing images.

I found a research article, "Learning Rotation-Invariant Convolutional Neural Networks for Object Detection in VHR Optical Remote Sensing Images", that is recognizing the objects in optical remote sensing images. They have used NWPU VHR-10 data set and its color images.

Therefore I started to check about the Optical Remote sensing images and I found this link. But in that link they didn't say anything about the color images or they have some different name (Multispectral imaging system, Panchromatic imaging system, and Superspectral Imaging Systems).

The color images used in this paper are Multispectral imaging system OR Panchromatic imaging system OR Superspectral Imaging Systems?
Sample image:


texts: Reducing and exporting to CSV returns duplicated rows with 0 values for calculated NDVI

texts:
I have table of 10 points that I want to extract the mean and standard deviation of NDVI from sentinel within a buffered area for a time period. I am running into a problem where when I reduce it and try to export it to a csv, I get a lot of duplicated rows with 0 values for the calculated statistics.

Here is the code I am working with:
print(sample);

// buffering around points
var bufferPoly = function(feature) {
  return feature.buffer(2500);   // adjust for buffer size in meters
};

var zone = sample.map(bufferPoly);

print(zone);

// Importing sentinel and cropping around feature collection
var sentinel = ee.ImageCollection("COPERNICUS/S2_SR").filterDate('2020-10-01', '2020-10-05')
    .select([
      'B4','B8',
    ])
    .filterBounds(zone)
    .map(function(image){return image.clip(zone)});

var count = sentinel.size(); // seeing how many bands there are...shouldn't be too many
print('Bands of sentinel: ', count);

print(sentinel); // we can see what is actually here, it is a collection of 12 images each having two bands

// Calculating NDVI for each image

// NDVI Function
function ndvi(im) {
  return im.normalizedDifference(['B8', 'B4']);
}

// mapping the ndvi function over all images
var ndvi = sentinel.map(ndvi);

var count = ndvi.size(); // seeing how many bands there are
print('Bands of ndvi: ', count);

print(ndvi); //We can see that it is a collection of 12 images with one ndvi band

//Now to reduce the image collection in a summary

//now extracting to the polygon

var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

var scale = ndvi.first().projection().nominalScale();

var reduced = ndvi.map(function(image){
  return image.reduceRegions({
    collection:zone ,
    reducer:reducers,
    scale: scale
  });
});

var count = reduced.size(); // There are 12 bands here
print('Bands of reduced: ', count);

print(reduced);

var table = reduced.flatten();

var count = table.size(); // There are 120 bands here
print('Bands of table: ', count);

print(table);

Export.table.toDrive({
  collection: table,
  description:'Example',
  fileFormat: 'csv'
});

For some reason, this is giving me 12 duplicated rows per the original row in the table. Any idea what is happening here?
link to code: https://code.earthengine.google.com/?scriptPath=users%2Fddlawton%2FAUS_SOILS_TEST%3AExtraction_code


texts:
I have given two co-ordinate pairs [{long, Lat}, {Long, Lat}]
like this i have several co-ordinate pairs.I need to create bounding box with each co-ordinate pair.

Now i have one co-ordinate pair[{Long, Lat}]

I have to find with this one co-ordinate pair whether this intersects any of the bounding box.

Using Nodejs i can do this easily

var bbox = new BoundingBox({ minlat: 48.123, minlon: 16.23, maxlat: 49.012, maxlon: 16.367 })
var bbox2 = new BoundingBox({ lat: 48.5, lon: 16.5 })
console.log(bbox.intersects(bbox2)) // true

But i have to store those co-rdinate pairs in postgis as bounding box and find the intersecting pair alone.

I have gone through certain links but i could not get clear knowledge out of it.


texts: Print City Map in High Resolution with TileMill (OSM Data)

texts:
I'm rather beginner in Print Map and GIS. I'd like to print the map of my little city using the data in OpenStreetMap. I've got an environment(local) with OSM DataBase  and I'm using TileMill to create StyleSheet and export data.
My problem is the definition of export. As you know, in TileMill you can export data in SVG, PDF and XML, and you can define just the BBox. The Resolution seems hard coded and set to 72dpi. To increase the quality and the details of export image I have to reduce the BBox (since to have a higher zoom level).

Since I need to print a classic city map (the paper size is, more or less, 80cm X 50cm and the city area 2km X 1km), does anybody have an idea for print it with a good resolution?

If i set 300 dpi as resolution, I'll solve my issue, but maybe I can't.
As I said before, I have a local environment, but perhaps someone knows some other ways.

Any ideas?


texts: FME Desktop 2013 Rastermosaiker: Why do I get strange edge patterns?

texts:
I have just mosaiked around 200 ecw tiles to a single ecw file. On inspecting the result, I see that I fme has produced a strange pattern around the edge. The Pixel values vary from R-G-B: 1 -1 -1  to 2-1-1 to 0-1-1 to 0-0-1and similar combinations very reaching over 2.

Can anyone offer an explanation and a solution so that this doesn´t happen?


texts: Attribute joins - string and text fields

texts:
I am trying to perform an attribute join, and it's just not happening for me.  The Excel field (csv) is formatted as text.  This is because it is prefixed with a 0 which is lost if it's any other format (e.g. 0100).  The field I am trying to join to in the attributes table is "string".

I am using ArcGIS 10.3.1

In the join pane, the table I am wanting to join (excel) is visible in drop down 2, with the fields visible in drop down 3.  However, when I select the field that the layer will join (attributes table) (drop down 1) all my field options in drop down 3 disappear.

What am I doing wrong?


texts:
I would like to calculate a buffer of 5km around country polygons. However, I somehow always get a warning and an error. Here is a reproducible example code.

library(sp)
library(sf)
library(tmap)
library(rnaturalearth)
countries_sp <- ne_countries(scale = 10)

countries_sp <- countries_sp[(countries_sp@data$geounit %in% "France"), ]

proj4string(countries_sp) <- CRS("+init=epsg:4326")

countries_for_buff <- st_as_sf(countries_sp, crs ="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

countries_buff_5km <- st_buffer( countries_for_buff, 5000)

here the warning:

dist is assumed to be in decimal degrees (arc_degrees).
Warning message:
In st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle,  :
st_buffer does not correctly buffer longitude/latitude data

Then

tm_shape(countries_sp) +
  tm_polygons() +
tm_shape(countries_buff_5km) +
  tm_polygons()

gives the error

 Error in if (all(tiny) && show.messages) warning("all relative area sizes are below thres.poly",  :
   missing value where TRUE/FALSE needed

What am I doing wrong? I guess there is something about the projection that does not quite work out the way it should.


texts: QGIS Geometry Generator not working for some features at some scales

texts:
I have a set of layers, one of them is labeled using geometry generator. Thing is, at some scales (1:800 and up) and some features (like 1/30) the geometry generator evalues to null geometry, so the feature is not labeled. I "fixed" it temporarily evaluating the resulting geometry and using another tye of label if it is null, but its not acceptable for printing. The Geometry Generator expression is
with_variable(
    'g',
    difference(
        buffer($geometry,12),
        aggregate(
            layer:='ManzanasM',
            aggregate:='collect',
            expression:=$geometry)),
    if(buffer( @g,-5) is null,
        buffer(pole_of_inaccessibility($geometry,1),5),
        buffer( @g,-5)))

It makes a buffer of itself and then a difference between that buffer and the layer 'ManzanasM', so the resulting polygon is outside the block. Its a grose replacement of a "label outside polygon" feature that I know its now implemented, but I cannot install new QGis in my PC without a lot of bureaucracy.
The layer is part of an Atlas dependant layout, I dont know if that's affecting it, but I dont think so. Used the same expression in Symbology to visualize exactly what im labeling, and after I zoom out more than 1:800 the geometry dissapears (and is replaced by a circle, because that's my temporary solution).
In the image, its visible how the colored parts outside the blocks are the zones where i want my label to be. There are two geometries that are not correctly created in the lower left.

EDIT:
Here are some example files. They need those custom expressions loaded.
https://drive.google.com/file/d/1XFJVQZfrADX92tnEe0Lz-k3AC9EGDFGA/view?usp=sharing


texts: How can I add a symbol to a callout label in QGIS?

texts:
I want to achieve something like this Barnett Bridge callout with a National Park Symbol underneath it.
I realize these maps are done using Illustrator, but I'm trying to stay in the realm of QGIS as much as possible.


texts: Estimating kernel bandwidth function with cross-validation gwrr package?

texts:
I am relatively new to R. I am attempting to use the gwrr package because I suspect that local collinearity may be an issue in my geographic weighted regression model.

If I am not mistaken, I first estimate the kernel bandwidth function using cross-validation. Below please find my R script. Diab is a spatial data set projected in Albers Equal Area Conic.

tt<-gwr.bw.est(dia2013~pctblacks + pcthis + pctpov+ lcollege+ newden+hosden10 +optden13 + totmdden13 + phyden13 + den13  +  nurseden13 + nohealthin + unemploy + transport, data=diab, kernel = "gauss", cv.tol=30)

However, I get the following error:

Error in (function (classes, fdef, mtable)  :

unable to find an inherited method for function ‘geometry’ for signature ‘"numeric"’

Does anyone know what is going on?


texts:
I am having trouble replacing the data source on CAD files using arcpy.

Everything I have tried runs into this error:

Runtime error type 'exceptions.ValueError' : Layer: Unexpected error

Here is a list of arguments I have tried:

layer.replaceDataSource('C:/CAD_DIR/CAD_File_1.dxf', 'CAD_WORKSPACE', 'Polyline')

layer.replaceDataSource('C:/CAD_DIR/CAD_File_1', 'CAD_WORKSPACE', 'Polyline')

layer.replaceDataSource('C:/CAD_DIR/', 'CAD_WORKSPACE', 'CAD_File_1/Polyline')

layer.replaceDataSource('C:/CAD_DIR/', 'CAD_WORKSPACE', 'CAD_File_1.dxf/Polyline')

layer.replaceDataSource('C:/CAD_DIR/', 'CAD_WORKSPACE', 'CAD_File_1.dxf.Polyline')

layer.replaceDataSource('C:/CAD_DIR/', 'CAD_WORKSPACE', 'CAD_File_1.Polyline')

It appears to take more time to think about going wrong when I put the 'CAD_File_1*' part in the 'Dataset_Name' section, that is to the right of 'CAD_WORKSPACE'.  But it still fails.

I have noticed that it treats the CAD file as a Feature Dataset and the particular feature type, in this case Polyline, as the Feature Class which may be the issue.

I have also tried this when pointing at a 'Polyline' layer:

layer.replaceDataSource('C:/CAD_DIR/', 'CAD_WORKSPACE', 'CAD_File_1.dxf)

It doesn't produce the same error or in fact any error but the data does not update.


texts: filter / generalize elevation points on the basis of some conditions

texts:
I wanted to generalize many elevation points in a rectangular area (about 700 square km).  I need to take about 3 points for square km. These points must satisfy some conditions:
One of these conditions is to find the highest point in a buffer of 500m: I would loop on all the points / "raster cell" (I think vector analysis will be too heavy, so I thought to do a "raster approach" creating an elevation grid from points) and give a score of "5" to the highest points / "raster cell" 500m around each.
I wanted to do it in an ArcGIS ModelBuilder flow.


texts: What is the easiest way to test python Add-Ins with ArcGIS 10.1 +?

texts:
I am trying to develop a Python Add-In for ArcGIS 10.3 and finding debugging pretty laborious. At the moment when I make a change to my python script I have to do the following:

Close ArcMap
Delete the "Python_Addins.esriaddin" file in my Default installation directory
double-click the "makeaddin.py" file to create another "Python_Addins.esriaddin" file
Install the add-in
Restart ArcMap and start the project
Test the tool

Surely, this can´t be the only way to test a python Add-In in ArcMap.


texts: How to convert Chinese Beijing 1954 coordinates to WGS84

texts:
I have a series of coordinates for locations in Sichuan Province China that I need to convert so that I can use them in Google Earth.  The coordinates are "1954 Beijing coordinates system" and show up as 3470081.60  36512300.88  with NO additional information other than elevation figures.

I ma pretty sure that the coordinate above corresponds with LAT 31.353305   LONG 108.131089 in WGS84 (Eastern Sichuan Province).

I'm struggling to understand what format the coordinates themselves are in - is it some type of modified MGRS?  I'm assuming there are two steps I need to follow here, firstly to convert the coordinate system itself and then secondly to covert from the Chinese Datum to WGS84....right?


texts:
Running ArcMap 10.3

I am displaying a layer of manholes from a table coming from SQL. The XYs are in decimal degrees however they are displaying no where near my other layers (NYC PLUTO). I have projected and reprojected the file but I want them to line up coming straight from SQL (not turning them into projected shapefiles). I made sure the data was correct as in the integers are correct (IE x -73.97 x  40.79).

I even created a new point layer, and calculated their centroids to make sure the numbers are similar (they are).

They are all in the same NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet projected coordinate system yet they refuse to go anywhere near each other.

I have a tons of data in SQL in the same dec degree format and would like to be able to easily display them all with parcel and street data.

Any insight on why these layers aren't lining up would be great.

UPDATE: using a geographic projection when displaying the points allowed SOME alignment however they are still off (manholes are not lined up on streets)


texts:
I'm fighting with an error I can't figure out where it is coming from.
Tried to compile all the things and used the latest unstable/stable ppa repositories as well, can't find the issue where it is comeing from.

Servers Setup is:

Postgres 9.1
Postgis 2.0 (2.1)
Mapnik 2.0
TileCache latest

An error occurred: Could not create datasource. No plugin found for type ' postgis
' (searched in: /usr/local/lib/mapnik/input) (encountered during parsing of layer 'landuse_gen0' in map '/var/mapnik/osm.xml')
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Service.py", line 343, in cgiHandler
    format, image = service.dispatchRequest( params, path_info, req_method, host )
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Service.py", line 208, in dispatchRequest
    return self.renderTile(tile, params.has_key('FORCE'))
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Service.py", line 138, in renderTile
    data = layer.render(tile, force=force)
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Layer.py", line 436, in render
    image = self.renderMetaTile(metatile, tile)
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Layer.py", line 396, in renderMetaTile
    data = self.renderTile(metatile)
  File "/usr/local/lib/python2.7/dist-packages/TileCache-2.11-py2.7.egg/TileCache/Layers/Mapnik.py", line 39, in renderTile
    mapnik.load_map(m,self.mapfile)


texts: Is it possible to perform automatic snapping of polygons in QGIS or ArcGIS?

texts: Note: in the picture below white space between the polygons is blank space – it does not contain any polygons as a result of Intersect or Union functions.


texts: I have a point layer for the cold air flow in specific region. I would like to simplify all the single points to a polygon shape as to be able to join the data with other spatial data. My way is pretty complicated (buffer the points e.g.400m, dissolve, simplify) because I have to extract the single values (orange, blue, light blue) and then do it with all these single values separately. Afterwards there is the need to create differences in order to avoid doubles. I wonder if somebody knows a more convenient way.


texts: QGIS: Python numpy error when installing plugin

texts:
I want to install the plugin "UMEP", but it results in this error:
An error has occurred while executing Python code:

RuntimeError: module compiled against API version 0xe but this version of numpy is 0xd
RuntimeError: module compiled against API version 0xe but this version of numpy is 0xd

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.16.4-Hannover Hannover, 654e76b3cb

Seems like a problem with numpy, I tried pip install numpy --upgrade in the windows command window, but it didn't solve the problem. You guys got any ideas?


texts: Troubleshooting File Geodatabase Parsing in Python / Fiona?

texts:
I have a test GDB file-folder that I have been able to verify that it does contain features through using some of ESRI's tools.

The code I have seems simple enough, but it keeps crashing when I try to call the next() function to get the first feature.  If I comment that function call out and just simply look through the gdb file for features, the program doesn't crash but prints no features.  What am I doing wrong?  Below is my code snippet:

gdb = fiona.open(GDB_Path)
gdb.driver # OpenFileGDB
gdb.schema # {'geometry': '3D MultiLineString', 'properties': OrderedDict([(u'FCSubtype', 'int'), ...
gdb.crs # {'init': u'epsg:4326'}

print "Getting first feature..."
gdb.next()

print "going to print all features..."
for feature in gdb:
   print feature
print "finished parsing GDB file"

Traceback shows that line 318 on fiona's collection.py file is crashing within the next call.  Other than trying another GDB file, is there anything else I should consider doing?


texts:
I have a huge amount of data in the form of CSV containing UTM coordinates as X and Y and an elevation value as Z information. I need to convert these data into a DEM as GeoTIFF for further analysis. In this case, a huge amount means 16 m. lines, with one point in X, Y and Z per line. The points are equally distributed, therefore an interpolation is not needed; each point just needs to be converted into a raster cell.

The original data came without separator, with fixed column widths. I already figured out how to convert the file syntax to use a separator instead of fixed widths and eliminate all space characters, using the stream text editor sed. From here on, normally my workflow would be to import the data into ArcGIS by creating a feature class from the X, Y and Z data and as a second step, converting the point shapefile into a GeoTIFF, using the Point to Raster tool. However, the file I currently have is way too big for this process.

Instead of the above described workflow, I was looking for an efficient alternative and discovered GDAL. However, in gdal_translate, the closest supported format I can find in the supported filetype list, is ASCII grid but no comma separated XYZ. Another difficulty is, that I have UTM coordinates, while most examples seem to use decimal degree coordinates. However, I need to stay within the UTM system (or at least, my output GeoTIFF needs to be in a UTM coordinate system).

So I am looking for a way to convert the CSV XYZ into a GeoTIFF, using GDAL, but so far wasn't able to find examples dealing with this exact problem. I would be very happy for some hints or even code examples.


texts: I've installed OpenGeo Suite, how do I use it with Leaflet?

texts:
I've just finished installing OpenGeo Suite on my VM Ubuntu web server.  I'd like to be able to use Leaflet from leafletjs.com to view my maps with on my mobile devices.  I'm supposed to create an HTML page, I assume, and use different snippets of code to make my map with it, however I don't know how to link the OpenGeo and Leaflet together.

It mentions using a tile layer URL, however I cannot figure out how to create this URL or where to find it. Also, do I have to create a new Leaflet HTML page for every map I want or am I way off? I'm sure I am missing some key elements of exactly how OpenGeo works, but I'm having a hard time finding any beginners guides and the documentation isn't very beginner friendly in some areas.


texts: Support materials and shaders in ArcGIS when importing model from 3D Max?

texts:
Does ArcGIS support materials and shaders when importing a model from 3D Max? Can I change the texture of the model when it is already added to the scene? How to do it?

Esri::ArcGISRuntime::SceneQuickView* m_sceneView = nullptr;
Esri::ArcGISRuntime::AnalysisOverlay* m_analysisOverlay = nullptr;
Esri::ArcGISRuntime::GeoElementViewshed* m_viewshed = nullptr;
Esri::ArcGISRuntime::GraphicsOverlay* m_graphicsOverlay = nullptr;
Esri::ArcGISRuntime::Graphic* m_tank = nullptr;

My addition of a 3D model to the scene:

void ViewshedGeoElement::createGraphic()
{
    // Create the Graphic Point
    const double x = -4.508708007847015;
    const double y = 48.38823243446344;
    const double z = 0;
    const Point tankPoint(x, y, z, SpatialReference(4326));
    const QString dataPath = QQmlProperty::read(this, "dataPath").toString();
    const float scale = 1;

    // Create the Graphic Symbol
    ModelSceneSymbol* sceneSymbol = new ModelSceneSymbol(QUrl(dataPath + "/2C19Msta-C/2c19Msta.3DS"), scale, this);
    sceneSymbol->setAnchorPosition(SceneSymbolAnchorPosition::Bottom);
    sceneSymbol->setHeading(180.0f);

    // Create the Graphic
    QVariantMap attr;
    attr[m_headingAttr] = 150.0;
    m_tank = new Graphic(tankPoint, attr, sceneSymbol, this);
    m_graphicsOverlay->graphics()->append(m_tank);
}


texts: Fixing invalid geometry (self intersection)

texts:
Attempted to extract raster by mask layer in QGIS 2.18.15 and for one of my shapefiles, the command could not be executed. The error message reads "Warning 1: Ring Self-intersection at or near point 97.520362854004077 18.494590759277401
ERROR 1: Cutline polygon is invalid."

I executed the Check Validity function, and identified the source of this self-intersection.

In trying to correct this 1 error, I have tried the following

(1) added a zero distance buffer using QGIS default plugin
(2) added a zero distance buffer using GRASS v.buffer.distance
(3) used multiple permutations of GRASS's v.clean & v.clean.advanced

But none of the above has rectified this issue i.e. when I run the Check Validity function, the error still persists. I have browsed through several threads within this forum but their solutions don't seem to work for me.