texts: Doing pairwise or pixel wise raster comparison using PyQGIS & Processing Scripts?

texts:
Do you know how to do a pairwise o pixel wise raster comparison using PyQGIS or Processing Scripts?

Pairwise or pixel wise comparisons are made between two raster datasets pixel by pixel. That means a pixel in an image a at a[n,n] will be compared with an image b at pixel b[n,n].

Comparisons like this are performed with functions like maximun(a,b)->value, minimun(a,b)->value, equal(a,b)->boolean, etc.

I have found and article that said by using the QgsRasterCalculator I can use functions like this:

le([rast1]@1, [rast2]@1, [rast2]@1)
or
gt([rast1]@1, [rast2]@1, [rast2]@1)

but none of these worked.

Are there any QgsRasterCalculator function I can use for performing pairwise analysis?

or

Do I have to convert my QgsRater to a Numpy matrix and to a QgsRaster again to perform this kind of analysis? (at first look seems a no so very efficient and CPU time consuming method)


texts: Problem running mvn install for geoexplorer

texts:
What's the solution to this: http://www.paste.to/MjU5ODE=

I'm getting the same error on Ubuntu 12.04.
Installed with: Java 1.6.0_27 (OpenJDK), Maven 2.2.1 and Ant 1.8.2

Is the Ant version the problem since 1.8.1 is in the maven repository?

How do I resolve this? Sorry I have little experience with Maven/Ant.


texts: Avalanche runout analysis within the latest TauDEM-package does not work (Error 000732: Parameters are not vaild)

texts:
I am trying to calculate the runout of avalanches with the TauDEM-toolbox (Version 5.3.7). I am running ArcGIS 10.4 and Windows 10.
I have tried using the sample data provided, following their tutorial and and it worked. When I am using my data (which I converted from .asc into .tif format with QGIS beforehand) I get the following message: ExecuteError: Failed to execute: Parameters are not valid.

All data have the same coordinate system as well as the same projection.


texts: GEE: Divide satellite imagery pixels into smaller sized pixels, without projection being changed

texts:
I would like to divide the moderate resolution satellite imagery into smaller pixels. To clarify this, I took a Landsat 30x30m image. For further analysis, I first want to split all pixels into 3x3m pixels, exactly corresponding with the original pixels. The method ee.Image.reproject() seems to come closest, although the pixels are reprojected into a slight different projection and therefore the pixels do not match with the original Landsat pixels. How can I Create an image where all new smaller pixels exactly match the original image?

Here a piece of simplified code which shows the mismatch (zoomed in to the corner of the image, where the mismatch is most clearly seen):

The mismatch removes when '.clip' is done at the end, where the image is added to the map

// Extent for the images to clip
var Haurvig = ee.Geometry.Rectangle([8.1247, 55.927, 8.1625, 55.907]);

// Load an LS8 image with least clouds
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                     .filterDate('2018-01-01','2018-12-01')
                     .filterBounds(Haurvig)
                     .sort('CLOUD_COVER').first();
// rename the band names and clip the image
var image =  image.select(['B2','B3','B4','B5','B6','B7','B10','pixel_qa'])
          .rename(['blue','green','red','nir','swir1','swir2','tir','qa'])
          .clip(Haurvig);

// Reproject the image to 3x3 m
var reprojected = image.reproject(image.projection().crs(), null, 3);

// add to the map
var trueColor = {bands: ['red','green','blue'], min: 300, max: 3800};
Map.addLayer(reprojected, trueColor, 'Reprojected image');
Map.addLayer(image, trueColor, 'Original image');
Map.centerObject(ee.Geometry.Point([8.162408851233522,55.907053640896]), 20);


texts:
I have multiple instances of my application for my users. Each user can specify an area which I want to load in my own PostgreSQL server, that make it's possible to render maps with my own mapnik stylesheet and nik4.

The idea was to load the data via overpass-api.de and load it with osm2psql. Each region has a size about 150MB. To separate the data, I want to use the --prefix option. Then I thought I could download the region each month again, use --append instead of --create option in osm2pql and the data will be updated. But as I know now, that is not possible.

A possible solution is, that I drop the tables and load it again. Another could be osmosis. If I don't delete the downloaded file I could create a changeset that is possible to load with the append option. But I have to save the old download file. There is a write write-pgsql-change option but it seems not to work with this multiple instances of the osm2psql import.

Is there any other solution to solve this problem or maybe I'm on a complete wrong way - so I'm open for any suggestions.

PS: I'm not able to load the whole planet, due to hardware limitations. At the moment, there are only 3 users, so the maintaining is quite small at the moment. But in future the will be a bit more users, and I don't know which region they choose, so I'm not able to download some geofabrik extracts and update them daily. And it seems much faster to handle much smaller tables with smaller extracts.


texts: I would like to use DEM to extract a river network but I don't want to use the polylines. Is it possible for me to extract polygon (or raster) rivers from DEM? I'm using both ArcGIS 10 and R. If possible, how reliable would the river width be?


texts: Convert multi-layer DWG into several KML or geoJSON files (one per layer)

texts:
I have some big DWG files, with 5 to 50 layers, which I draw and geolocated using AutoCAD (with some external modules).
 I need to convert them to KML or geoJSON so that the JS library can load them and display them on the map.

I used AutoCAD to export my .DWG into .DXF, then I installed GDAL and used ogr2ogr to convert those DXFs into geoJSON. It works, but the generated geoJSON is very big (more than 10 megabytes) and it loads forever once I try to display it into a google map or a leaflet map.

What I would like to do, is to have a separated geoJSON/KML for every layer of my DWG/DXF file, so that I can load and display each layer separately which hopefully would increase performance.

I am trying to have an automated workflow, because as I have a lot of files, it would take too much time to take every one of them and split layers manually.

Any idea ?


texts: When are .prj files included with Shapefiles?

texts:
I notice that sometimes Shapefiles that I have received contain .prj files. Sometimes they don't. Wikipedia notes that the .prj file is not "mandatory" and that it contains the :

  projection format; the coordinate system and projection information, a plain text file describing the projection using well-known text format

So my question is: why is it sometimes included, and sometimes not? Do newer versions of Esri tools include it, whereas older ones didn't?


texts:
I have a raster layer with pixel values of 1 and 0. I would like to extract values from several other, overlaying, raster layers for all the cells that equal 1.

Is there any ArcMap tool similar to extract multi values to point that I could use?


texts:
This leaflet fiddle produce what I want. eg. tooltip on custom marker

https://playground-leaflet.rhcloud.com/juma/edit?html,output

I try to implement the same solution when I am fetching my data from an ajax call

I tried this, but it does not work, the console log is telling me Failed to execute 'appendChild' on 'Node': parameter 1 is not of type 'Node'.(…)

If i try to bindTooltip as text it works, but not when I try to reach the properties for each feature

    var myicon = L.icon({
                            iconUrl: 'assets/img/point.png',
                            iconSize: [8, 8]
                            });

    var pointlayer= L.geoJson(null,{
    pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {icon: sandicon});
    },
    onEachFeature: function (feature, layer) {
        layer.bindTooltip(feature.properties.gid, {permanent:true,sticky:true}).openTooltip();
    }
    }).addTo(map);

    $.ajax({
        url: "root",
        dataType: 'jsonp',
        cache: false,
        jsonpCallback: 'getgjson',
        success: handledata
    });

    function handledata(data) {
        pointlayer.addData(data);
    };


texts: Display more than one attribute in symbology in ArcGIS Pro?

texts:
I want to color code temperature and represent precipitation by dot density in the same map.

How can I do this in ArcGIS Pro?


texts: Splitting multiple noncontiguous polylines into multiple files

texts: I have multiple polyline shapefiles that are non-contiguous, that is to say that each shapefile contains polylines that have breaks in them. I am looking to split each of these polyline segments into its own shapefile. Furthermore, I need to do this for several hundred shapefiles, so I am wondering if there is a code I can write in python, or else a function that exists in ArcGIS itself to do this quickly for all of the shapefiles.


texts: QGIS Legend only with Raster Minimum and Maximum values

texts: I want to have a mimim and maximum raster values in the legend but it appears to me a long list of all values how can I resume the legend minimum and maximum values:


texts:
Where do I find the MapInfo Pro settings file?

I just opened program and it looks like after installation, all my changes are back to default.

I'm using MapInfo Pro 15.0 32-bit


texts: Is there any way to get the SQL View parameter list via REST API?

texts:
I have a view SQL layer with some parameters. I need to use the REST API to guess these parameters and show to the user in the interface to he/she fill the values before submit.

Going to rest/workspaces/<WS>/datastores/<DS>/featuretypes/ I can see the layer attributes but no info about the view SQL.


texts:
I am using the qgis python editor to create a script and run it. When I am done writing my script, I would like to save it. However when I click on Save as I get the following error

An error has occurred while executing Python code:

AttributeError: 'QDir' object has no attribute 'replace'
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python\console\console.py", line 660, in saveAsScriptFile
    pathFileName = os.path.join(folder, fileName)
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\ntpath.py", line 65, in join
    result_drive, result_path = splitdrive(path)
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\ntpath.py", line 116, in splitdrive
    normp = p.replace(altsep, sep)
AttributeError: 'QDir' object has no attribute 'replace'

Python version: 2.7.14 (v2.7.14:84471935ed, Sep 16 2017, 20:25:58) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.22 Las Palmas, 5a99678469

Python Path:
C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python/plugins\processing
C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python
C:/Users/aboufira/.qgis2/python
C:/Users/aboufira/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python/plugins
C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.18\bin
C:\PROGRA~1\QGIS2~1.18\apps\Python27
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/aboufira/.qgis2//python
C:\Users\aboufira\.qgis2\python\plugins\DigitizingTools\tools
C:/Users/aboufira/.qgis2/python/plugins\pdokservicesplugin

I have no idea how to resolve this, has anyone seen this before?


texts: Listing geodatabases in local SQL Express Workgroup ArcSDE using ArcPy/Python?

texts:
I would like to create a list of the multiple geodatabases that sit in my local SQL Express Workgroup SDE server.

I have set up a connection using the "Database Connections\Connection to MyComputer_SQLEXPRESS.sde", but there doesn't seem to be any arcpy functions that expose the geodatabases.

It seems to identify the connection to a workspace, but doesn't give any information for what is in the workspace.

I know I can connect to each geodatabase separately and retrieve information about each of them, but I really need to get just a list of the geodatabases themselves.


texts: canvasMoveEvent: canvas not refreshing in pyqgis

texts:
I wrote a qgis plugin where I can modify some features interactively. I have a class MovePointTool derived from QgsMapToolIdentify.

When a point is moved, I call a function that does some (relatively light) processing and modifies some other dependent layers (self.database.moveIP and self.database.display in the snippet below).

class MovePointTool(QgsMapToolIdentify):

    def __init__(self, mapCanvas, layer, database):
        self.ready = True
        print("Creating move instance")
        QgsMapToolIdentify.__init__(self, mapCanvas)
        self.setCursor(Qt.CrossCursor)
        self.layer    = layer
        self.dragging = False
        self.feature  = None
        self.database = database

    def canvasPressEvent(self, event):
        print("press")
        found_features = self.identify(event.x(), event.y(),
                                       [self.layer],
                                       self.TopDownAll)
        if len(found_features) > 0:
            self.dragging = True
            self.feature  = found_features[0].mFeature
        else:
            self.dragging = False
            self.feature  = None

    def canvasMoveEvent(self, event):
        if self.dragging and self.ready:
            self.ready = False
            point = self.toLayerCoordinates(self.layer,
                                            event.pos())
            geometry = QgsGeometry.fromPoint(point)
            self.layer.changeGeometry(self.feature.id(), geometry)
            self.database.moveIP(self.feature.id(), self.toLayerCoordinates(self.layer, event.pos()))
            self.database.display()
            self.canvas().refresh()
            print('c')
            self.ready = True

    def canvasReleaseEvent(self, event):
            print("release")
            self.dragging = False
            self.feature  = None

I am not sure how QgsMapToolIdentify works, and if it is asynchronous. I have no experience with asynchronous programming. I had to add the self.ready flag above to avoid segmentation faults (not sure how robust it is).

However, the canvas does not refresh smoothly (lag and miss some moves), except if the

print('c')

line is present and IF the python console if open. Then, it refreshes very smoothly.

This is a strange issue. What is the proper way to be able to move a point, do some processing based on this move and then refresh the canvas? Is there any barrier to add somewhere?


texts: Generating new shp files according to attribute in QGIS

texts: I know that it is possible to color a vector map by attribute. I separated a map into 4 colors based on an attribute and now I want to generate a new map for each color. Is there a way to do this in QGIS?


texts: Rendering order of layers in QGis libraries

texts:
I am using the QGis libraries in a external Qt App.
When adding vector and raster layers (in any order) to the QgsMapLayerRegistry instance, it somehow sorts the layers so that the vector layers is first in the QgsMapLayerRegistry.mMapLayers list and then the raster layers.
Which make perfectly sense otherwise at render time you wont see the vector layers if the order was the other way round.

But when adding a vector point layer with a memory data provider, i would like the layer to be first in the list, or at least "sorted" with the other vector layers.
But is always end up last in the QgsMapLayerRegistry.mMapLayers list, and at render time it always gets rendered "beneath" the raster layer, how can i add the vector point layer with the memory data model, to be rendered "on top" of all other layers?


texts: Plotting two shapefiles in same plot. Legend of 1st shapefile hidden by second shapefile. Python GeoPandas

texts:
With
ax = gdf.plot(cmap="jet", column=ColumnaAVisualizar, categorical=True,figsize=(12,12),legend=True)
ax2=ZDE.boundary.plot(ls="--",color="k", label="Zona de Estudio", ax=ax)
plt.show()

It gives:

Now with
ax = gdf.plot(cmap="jet", column=ColumnaAVisualizar, categorical=True,figsize=(12,12),legend=True)
l1=ax.legend(loc='center right', bbox_to_anchor=(1.0, 0.5))
ax2=ZDE.boundary.plot(ls="--",color="k", label="Zona de Estudio", ax=ax)
l2=ax2.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
ax.add_artist(l1)
plt.show()

It gives:

So I would like something with the two legends:

PS: Both gdf have the same CRS.


texts: Drawing median line from 2 points on Google Earth?

texts:
Given 2 locations A and B, I am looking to draw a line around the earth where for any point C on the line, length AC = length BC on the earth surface.
That is to say the line where we always are at the same distance from A and B.
I would like to be able to draw such a line on Google Earth.
I thought this was going to be easy but after an extensive research, I couldn't find any KML example to do that.


texts: Can anything be inferred if NHD doesn't name a stream but the local municipality does?

texts:
Case in point: https://watersgeo.epa.gov/watershedreport/?comid=5781291 It shows up in the NHD as not having a name but I'm pretty sure that the local municipality that that creek is in (City of Austin, Texas) calls it the Fort Branch Creek: http://austintexas.gov/sites/default/files/files/Watershed/eii/Fort_EII_ph1_2011.pdf
Maybe creeks need to be of a minimum length to be named in the NHD?


texts:
I have created a shapefile with approximately 300 feature of shapeType POINTZ. I can fill my attributs and define point.
Here is my function that i call for each feature:

def upd_shp(fn, geo, attr):
    r = shapefile.Reader(fn)
    w = shapefile.Writer(r.shapeType)
    w.fields = list(r.fields)
    w.records.extend(r.records())
    w._shapes.extend(r.shapes())
    x = geo[1]
    y = geo[0]
    z = geo[2]
    w.point(x, y, z)
    w.record(*attr)
    w.save(fn)

Then when i do ogrinfo -al myshp.shp i can see POINT (8.833808735269852 42.540690712189829) but there is not  'z' value.
When i do ogrinfo -al example.shp which it's a shapefile made by my company i can see POINT (8.75044460305986 42.567946874150579 37.384099928078562).

I know my 'z' is in the shapefile but i don't understand why i can't see it with ogrinfo.

I searched and I found Creating and reading z values (coordinates) with pyshp? but it did not help me much.

I would like to know if i'm missing something?


texts: ArcGIS API Dojo problem to activate a button for a function

texts:
I'm trying to connect (re-connect) a function to a panel button the dom by id is ok, if I use a function like :

on(selectPoly, "click", function(){
   html.set(dom.byId("drawToolTips"),"Autre message !");
})

It work well but when I try the code below I can't reach my start select function (the other function who are not in a panel work well)

var selectPoly = dom.byId("selectPoly");
var selectCircle = dom.byId("selectCircle");
var selectExtent = dom.byId("selectExtent");

var validation = dom.byId("end");
var eraseAll = dom.byId("del");
var goBck = dom.byId("goBack");

// Connect bootstrap tools buttons
on(selectPoly, "click", startSelect("Polygon"));
on(selectCircle, "click", startSelect("Circle"));
on(selectExtent, "click", startSelect("Extent"));
on(validation, "click", endSelect);
on(eraseAll, "click", delSelect);
on(goBck, "click", goBackfn);

function startSelect(drawType) {
  console.log(drawType);

  //activate module zoning séléction
  switch (drawType) {
    case "Polygon":
      select.activate(Draw.POLYGON);
      html.set(dom.byId("drawToolTips"),
        "Cliquez pour ajouter un point,
 \n\
            Double-Cliquez pour Terminer la forme,\n\

 Echap pour annuler le tracé en cours");

      // fallthrough
    case "Circle":
      html.set(dom.byId("drawToolTips"), "Autre message !");
      // select.activate(Draw.CIRCLE);
      break;

  }
};

here is my html panel :

        Séléction géographique

Cliquez sur un outil pour commencer la séléction


texts:
Background:
At 3.18.2 I have a Geopackage attribute table that I want to refactor, as shown in the screenshot below.  I want to limit the field box_id length to 8.
The table currently has 55 records, but will soon have tens of thousands.

Problem:
The refactoring correctly adjusts the field length.  However, it also changes the fid field from Int64 to Int8, as shown in the screenshot below.  Why is this happening, and how can I prevent it?  Neither the QGIS docs, nor SE GIS, address this problem.


texts: Entering Track from Garmin to GIS not getting basemap and orthomap to line up with Track?

texts:
I collected a GPS path on a flight and I would like to condense the path to the place where we landed and took pictures (condensed area at base of map; I want to delete the other part of the line). Ultimately, I would like to eliminate all of the information except for the location where the flight landed.

When I add the Track to GIS and "zoom to layer," it makes the data frame leader go off the page and the Track does not match up with the basemap.

What is the issue?


texts:
I would like to study a range of data using the same technique in ArcMAP (it will take too much time to do it one by one ;) !). I've been just recently playing with ArcMAP and Python (which can help me to generate an automatic sequence).

I also read that it is easier to import .csv file into ArcMAP.

I'm importing the csv file using csv.reader (just coordinate of points without header, two colums x and y) like:

file = csv.read(...)

Then I didn't managed to print the file, because i obtained this message:

<_csv.writer object at 0x1426F7B0>

And I want to put the first column in the x coordinate and the second column in the y coordinate in a point shapefile?

I read a lot of information how to create file. But I'm stuck in taking the first and second column of my file.

Can you give me some advices?

TA


texts: Installing LizardTech GeoExpress Command Line Applications on Ubuntu 16.04

texts:
I need to convert MrSid files so I can view them as TIFFs in QGIS.  I am trying to install the LizardTech GeoExpress Command Line Applications available here to do the conversions.  I have downloaded the correct Linux 64-bit version, and extracted the contents to a directory on the hard-drive.  The included read-me file states this...These utilities dynamically load libraries that are also installed.
Please ensure that your LD_LIBRARY_PATH variable includes the directory to
which this package is installed.

I am new to Linux and found this command online to map to the tools:

export LD_LIBRARY_PATH="/media/gerry/02689296-24e3-4e9b-b940-b33bdbbcab43/GIS/OtherSoftware/GeoExpressCLUtils-9.5.0.4326-linux64"

When I try to run a simple example I get this...

mrsidgeoinfo ortho_1-1_1n_s_wa073_2013_1.sid
mrsidgeoinfo: command not found

I am not too clear on on how to install these utilities. How do I install the LizardTech GeoExpress Command Line Applications?


texts: ArcGIS Online "Add members using their enterprise ID" - what value to use for "Enterprise ID"?

texts:
We have ArcGIS Online configured for "enterprise" logins to authenticate against our MS Active Directory.  I can add users from AD by sending them invitations in ArcGIS Online.  They must respond to these emailed invitations for their account to be created in ArcGIS Online.

I would like to add AD users using the option, "Add members using their enterprise ID without sending invitations".  This gets their account set up immediately, and then I can add them to groups, etc, without waiting for them to respond to an invitation.

However, to use this option, I have to include the user's "Enterprise ID".

What is this Enterprise ID, and how do I find out what a particular user's Enterprise ID is?  Is it something I have to get from Active Directory?  If so, how do I find it for various users?

(I have tried various combinations of username, email address, and organisation name and underscores, but with no luck.)

UPDATE

As per comments, I have now tried populating the "Enterprise ID" field using the format, DOMAIN\user.name (and DOMAIN/user.name), but it will not accept anything with either a \ or a / in it (no error, but the "Next" button does nothing).  As soon as I remove the (back)slash, the "Next" button works as expected (but the account does not work because I don't know what I'm supposed to put in the "Enterprise ID" field.

I note that the examples in ESRI's documentation do not use this format.  The examples there are "jonc1111" and "sati3554".  See:  https://doc.arcgis.com/en/arcgis-online/administer/invite-users.htm#ESRI_SECTION1_D4222EB63EF14C96BF37F1A098F37AFA

I feel like there's something very simple/obvious that I'm missing, like I'm looking at the problem through the wrong lens or something.  I just can't figure out what I'm sure is supposed to be straightforward.


texts: Classifying DEM according to elevation in R?

texts:
r<-raster("DEM.tif")

How to do the classification based on every ten meters for example and have different colors too?


texts: Gdal version 3s versus 2s: awakward result

texts:
I have a polygon in a kml file.
I wanted to create bounding box of this polygon and save as kml file.
The following code used to work with gdal version 2s.
But with gdal version 3s, the output kml file shows ocean!
The input kml is certainly belongs to land, far away from ocean.

import os, glob, ogr, osr, numpy as np
from shapely.wkt import loads
from shapely.geometry import Polygon

fi = 'AOI.kml'
ds = ogr.Open(fi,0)
slayer = ds.GetLayer(0)

fo = 'AOI_bound.kml'

ods =  ogr.GetDriverByName('kml').CreateDataSource(fo)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

olay = ods.CreateLayer('poly', srs, geom_type=ogr.wkbPolygon)
olay.CreateField(ogr.FieldDefn('Name', ogr.OFTString))

for feat in slayer:
    poly = loads(str(feat.GetGeometryRef())).envelope

    ogeom = ogr.CreateGeometryFromWkb(poly.wkb)
    ofeat = ogr.Feature(olay.GetLayerDefn())
    ofeat.SetGeometry(ogeom)
    ofeat.SetField('Name', 'poly')
    olay.CreateFeature(ofeat)
ds = ods = None


texts: QField Capture screen write text in a field

texts: Just started using QField with Samsung tablet. The tablet has a "screen write" function. Is there a way of linking this text to points in a QGIS/QField table? It would be more flexible than typing into a text box.


texts: Installing ExifRead module into QGIS 2.14

texts: I'm attempting to install ExifRead into QGIS. Assuming zero knowledge of python coding. Is there a command code I can send which will install the .whl? I'm having zero luck following python's own guide.


texts: Orfeo Toolbox (OTB) has a tool that I'd like to use for some image processing as part of a longer Python workflow that I'm executing in Google Colab. While OTB can be called via Python if the whole program is installed locally, I can't figure out if there is a way to install OTB's Python tools in a Colab notebook.


texts: Force virtual layer to based on layer's filter (query builder)

texts:
I have a layers, which are filtered by query builder.

Next I want to create virtual layer (VL), that is based on that filtered layers.
Sadly, I'm getting VL with all features from the layers. Query builder is ignored.

Is there a way to force VL to use an output from query builder, and not all features?


texts: Using nearest neighbor query in SpatiaLite?

texts:
I am trying to find the nearest line to each point in a layer, and I have over 1 million points. I want to be able to associated the closest unique line ID with the point unique IDs.

I have tried to use the solution given in this post Finding nearest neighbors between two tables with point locations in SpatiaLite? on a small subset, and I do get a result table but the distance value is null, and the IDs from the lines that are associated with each point are not the IDs for the closest lines.

I am new with SpatialLite and running spatial queries in SQL. I don't really understand the purpose of this statement: WHERE g1.OGC_FID <> g2.OGC_FID.

I have also played around with this line: AND ST_Contains(ST_Expand(g1.geometry,50),g2.geometry) as well as removed it and still get no distance values, even though I am getting an ID.

I was wondering if it might be a coordinate system issue but both the points and the lines are in the same system. I am using the QSpatiaLite plugin in QGIS. Can anyone provide some guidance?

Using the answer below I was able to get the query working on some test data. However I am having trouble applying the spatial index. I've tried the query below, but it runs slower than the query statement without the spatial index. What am I doing wrong?

UPDATE UgandaBikeTaxis_test SET dist=(SELECT
    ST_Distance(p.geometry, l.geometry) AND
p.ROWID IN
(SELECT ROWID FROM SpatialIndex WHERE f_table_name= "ugandabiketaxis_test"  AND search_frame=l.Geometry)
    FROM osm_filtered_prj_test AS l, ugandabiketaxis_test AS p
    WHERE p.uniqueID = ugandabiketaxis_test.uniqueID
    ORDER BY ST_Distance(p.geometry, l.geometry) LIMIT 1);


texts:
I'm trying to clip a set of polylines to a polygon with geopandas but it's returning an empty set. Here's my code:
>>> import geopandas as gpd
>>> import matplotlib.pyplot as plt

>>> gdf_lines = gpd.read_file('by_christina_nhd_0.geojson')
>>> gdf_basin = gpd.read_file('christina_proj_4269.geojson')

>>> ax = gdf_lines.plot()
>>> gdf_basin.plot(ax=ax, color='white', linewidth=2, edgecolor='black')
>>> plt.show()

>>> gdf_lines.crs == gdf_basin.crs
True

>>>clipped_lines = gpd.clip(gdf_lines, gdf_basin)
>>>clipped_lines.empty
True

They obviously intersect, so what's up? Why isn't clip returning those lines within the polygon? Any ideas?
Here is a link to the geojson files
FYI - I've also tried this on QGIS and I get an empty set there too. It also has a warning about "invalid geometry." Not sure if that has anything to do with it.


texts:
I have a CSV file with number of records that have WKT Point geometry.  I am bit lost on how to insert to this into a new point layer.

Sample data from CSV file from GeometryItem Field:

POINT (-13174097.3521 4006999.38409)
POINT (-13172281.9726 4006836.9554)
POINT (-13172215.0902 4009803.66756)
POINT (-13173987.4738 4011721.75923)

arcpy.env.overwriteOutput

out_path = r'D:\My Documents\Metro_RTSO\GISData\Metro_RTSO.gdb'
arcpy.env.workspace = out_path
sr = arcpy.SpatialReference(3395)
print sr.name
newFC= 'ITS_20171024_CCTVs'

if arcpy.Exists(newFC):
    print newFC + " Exists, Deleting"
    arcpy.Delete_management(newFC)
    arcpy.CreateFeatureclass_management(out_path, newFC, 'POINT', spatial_reference=sr)

else:
     arcpy.CreateFeatureclass_management(out_path, newFC, 'POINT', spatial_reference=sr)

with open(r'D:\My Documents\Metro_RTSO\GISData\ITS Inventory\ITS_20171024_CCTVs.csv', 'r') as f:
     reader = csv.DictReader(f)
     insertCursor = arcpy.da.InsertCursor(newFC, ['SHAPE@WKT'])
     Lines = 0
     records = 0
     for row in reader:
         Lines += 1
         if row["GeometryItem"] != None:
            records +=1
            ##print row["GeometryItem"]
            print arcpy.FromWKT(row["GeometryItem"],sr)
            new_pt = arcpy.Point(arcpy.FromWKT(row["GeometryItem"]))
            insertCursor.insertRow(new_pt)
     del insertCursor
     print records, Lines

print 'finished'


texts: Feature Grid Layer to cluster features at higher zoom levels in Leaflet

texts:
I was wondering if there is a plug in for Leaflet which accomplishes clustering features within a grid layer. Similar to marker clusters this allows for rendering an abstraction of all features with a count of the features within that area. In my case I want to do this for centerlines of a city which is too much to render at higher zoom levels.

The overall end effect is similar to what is seen in the LA geohub:

where you need to click the grid till at proper zoom level to render any features.


texts: I managed to crash my ArcCatalog trying to bring multiple feature classes over into a new file geodatabase \ feature dataset. I tried starting up the program again and during the start up it crashed. I have restarted my computer yet ArcCatalog will still crash upon start up. I am unsure how to solve this problem apart from uninstalling ArcGIS for Desktop which I do not want to do. Has anyone had any experience with this issue and/or know a way to only re-install ArcCatalog?


texts: Is it possible to increase Output Field Name length beyond 10 characters?

texts:
When creating a new Column, the default value for the length of the name is 10. Using the Field Calculator via the Model Builder, I want to include the names of regional areas into the columns which has more than 10 characters.

Can the length of the Output Field Name be increased?


texts: 64 bit ArcPy works on server but not through Citrix, path looks fine, what is going on? (ArcGIS Desktop 10.3)

texts:
I can import ArcPy fine if I'm on the server itself; but If I access ArcGIS through Citrix using the same server I can't import it. I checked the paths and they look identical, I don't see anything missing.

I noticed that on the actual server, there is no Python path; while on the Citrix applicaiton, the path is pointing to regular 32 bit

>>> print os.environ['PYTHONPATH']

c:\program files (x86)\arcgis\desktop10.3\bin;c:\program files (x86)\arcgis\desktop10.3\ArcPy;c:\program files (x86)\arcgis\desktop10.3\ArcToolbox\Scripts

This area is not my area of expertise, but I am guessing that the Python path needs to change for the Citrix application?


texts:
I have a line dataset that has road names and different surface types along different segments of the road. I want to calculate the total length of each surface type for each road.

I've been able to get the groups and counts as I need using virtual layers / SQL

Select Road, SURF_TYPE, count(SURF_TYPE), geometry from road_extract group by Road,SURF_TYPE

Get's me

Rather than geometry I want to get the total length of each surface type.

I've found one workaround where if I have a column with $length in the the source table I can include sum(length) in the SQL to get the total length for each surface type but is it possible without having this column?

Also is it possible to get this table formatted with the just the lengths like...

Road,    Sealed, Unsealed, Floodway, Total
Airport, 15    , 0       , 0       , 15
Balonne, 20    , 5       , 0       , 25
Burke,   40    , 20      , 2       , 62

I have tried the following with no luck

select
 Road, SURF_TYPE,
 sum(length)
from
 road_extract

 Group by
  Grouping sets (
   (Road, SURF_TYPE, sum(length)),
   (Road),
   (SURF_TYPE),
   ()
  );


texts:
Does anyone know if there's a tool in QGIS that is equivalent to ArcGIS Point Statistics?

I need to create a raster with values from an attribute in a Point layer.


texts: Geoserver SLD - priority option displays multiple labels

texts:
I'm using priority option for labels in my geoserver SLD stylesheet because I have multiple overlapping points.

<Priority>
  <ogc:PropertyName>language</ogc:PropertyName>
</Priority>

The problem is, geoserver renders both labels if there are two overlapping points. The priority one is above, and the another one is below.

I need only one of those labels (top priority one) to be rendered.

EDIT:

Whole TextSymbolizer looks like this:

<TextSymbolizer>
        <Label>
          <ogc:PropertyName>pisanjeimena</ogc:PropertyName>
        </Label>
        <Priority>
          <ogc:PropertyName>language</ogc:PropertyName>
        </Priority>
        <VendorOption name="group">false</VendorOption>
        <VendorOption name="labelAllGroup">true</VendorOption>
        <VendorOption name="spaceAround">10</VendorOption> <!--allow labels to overlap-->
        <!--<VendorOption name="followLine">true</VendorOption>--> <!--option forces a label to follow the curve of the line-->
        <VendorOption name="maxDisplacement">30</VendorOption> <!--option controls the displacement of the label along a line-->
        <VendorOption name="repeat">100</VendorOption> <!--option determines how often GeoServer displays labels along a line-->
        <VendorOption name="autoWrap">150</VendorOption> <!--option wraps labels when they exceed the given width (in pixels)-->
        <!--<VendorOption name="conflictResolution">false</VendorOption>--> <!--conflictResolution option to false causes this label to bypass conflict resolution-->
        <VendorOption name="goodnessOfFit">0.3</VendorOption>
        <Halo>
          <Radius>2</Radius>
          <Fill>
            <CssParameter name="fill">#FFFFFF</CssParameter>
            <CssParameter name="fill-opacity">0.6</CssParameter>
          </Fill>
        </Halo>
        <Font>
          <CssParameter name="font-family">Arial</CssParameter>
          <CssParameter name="font-size">12</CssParameter>
          <CssParameter name="font-style">normal</CssParameter>
          <CssParameter name="font-weight">bold</CssParameter>
          <CssParameter name="fill-opacity">0.6</CssParameter>
        </Font>
        <LabelPlacement>
          <PointPlacement>
            <AnchorPoint>
              <AnchorPointX>0.5</AnchorPointX>
              <AnchorPointY>0.0</AnchorPointY>
            </AnchorPoint>
            <Displacement>
              <DisplacementX>0</DisplacementX>
              <DisplacementY>5</DisplacementY>
            </Displacement>
          </PointPlacement>
        </LabelPlacement>
        <Fill>
          <CssParameter name="fill"><ogc:PropertyName>boja</ogc:PropertyName></CssParameter>
        </Fill>
      </TextSymbolizer>

Language property that is used in <Priority> is numerical value (language ID). Language with higher ID should be a priority one.


texts: PostGIS view on multiple geometry columns

texts:
I have two PostGIS tables, ‘filelinks’ and ‘maricopa_plss_qtrqtrs’. The filelinks table has a geometry type column named ‘geom_p’ and the maricopa_plss_qtrqtrs has a geometry type column named ‘geom’. The filelinks table can contain its own polygon geometry in geom_p, but if that geom_p value is null, then its ‘parent_ufid’ value will contain the ‘ufid’ value from maricopa_plss_qtrqtrs, which always has a populated geom field.

What I want to do is create a view (‘recursive_view’) that pulls the geometry from either the filelinks table geom_p (if it’s not null) or the maricopa_plss_qtrqtrs table’s geom value, related through the parent_ufid. Here’s what I have as the current view definition:

CREATE OR REPLACE VIEW gis.recursive_view AS
SELECT
    f.ufid,
    qq.geom,
    f.geom_p,
    f.name,
    f.uri,
    f.type,
    f.date
FROM gis.filelinks f, gis.maricopa_plss_qtrqtrs qq
WHERE qq.ufid = f.parent_ufid or f.geom_p is not null;

The error I’m getting is:

ERROR:  cannot change name of view column "geom_p" to "geom"
********** Error **********
ERROR: cannot change name of view column "geom_p" to "geom"
SQL state: 42P16

I thought a PostGIS view can have multiple geometry columns, just like a table, no? What is the problem with my view definition then?


texts: Creating geodatabase administrator account other than SDE?

texts:
I need to create another account other than SDE and make it geodatabase administrator. Is it possible?

I have already created another account and have given SUPERUSER priviledge. However I cannot do GDB admin tasks like database compression.

Do I have to change ArcSDE system tables too?

I am working with PostreSQL RDBMS.


texts:
I'm very new to postgis and pgRouting but I'm using it for my thesys.
I need any map to train my self with routing and TSP problems in pgRouting, but on the Internet I found unworking links only!
Can anyone suggest me where to download these maps (rather in osm)?
thanks a lot
Laura


texts:
I want to use OGR in Python to write a simple point-in-polygon test:

def point_in_polygon(point, polygon):
    """
    point : [longitude, latitude]

    polygon : [(lon1, lat1), (lon2, lat2), ..., (lonn, latn)]

    """
    # Create spatialReference
    spatialReference = osgeo.osr.SpatialReference()
    spatialReference.SetWellKnownGeogCS("WGS84")
    # Create ring
    ring = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
    # Add points
    for lon, lat in polygon:
        ring.AddPoint(lon, lat)
    # Create polygon
    poly = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
    poly.AssignSpatialReference(spatialReference)
    poly.AddGeometry(ring)
    # Create point
    pt = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
    pt.AssignSpatialReference(spatialReference)
    pt.SetPoint(0, point[0], point[1])
    # test
    return pt.Within(poly)

However, it doesn't seem to work properly:

In [22]: polygon = [(-10., -10.), (-10., 10.), (10., 10.), (10., -10.), (-10., -10.)]

In [23]: point_in_polygon([0., 0.], polygon)
Out[23]: True

In [24]: point_in_polygon([359., 0.], polygon)
Out[24]: False

Any ideas what I'm missing?


texts: What's the difference between hypsography and topography?

texts:
I've recently come across some data that has been labeled "hypsography". I looked it up and it sounds like topography. But this data source has data labeled "topography" separate from its "hypsography". What's the difference?

I've read the wiki for 'Elevation' and its starting to sound like all the "topos" I've been reading are really "hypsos"


texts: SHP file Left / Right Label Alignment issue in QGIS 2.0.1-Dufour

texts: cld anyone let me know how it works to get the desired label placement for the opened SHP file into QGIS 2.0.1- Dufo


texts: R stars package and set dimensions to a NetCDF file

texts:
Related: GDAL and reprojecting a NetCDF file

Following above, I would like to reproject a NetCDF ( NOAA GOES-16 file ) using R stars package.
However, I have been unsuccessful at the first step; I mean, unable to assign original CRS information to the file.

Parameters I want to assign were obtained by gdalinfo:

Proj.4 string = "+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"
Origin = (-3627271.328111211769283,1583173.512170388596132)
Pixel Size = (1002.008657743770527,1002.008657743770527)

My failed R code is:

library(stars)              # stars package 0.3-0

nc1= "./OR_ABI-L1b-RadC-M3C01_G16_s20190611802132_e20190611804505_c20190611804549.nc"
(nc1_s1= read_stars(nc1, along= c("DQF", "Rad")))
crs_proj4 = "+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"

nc1_out = st_set_dimensions(nc1_s1,
                            offset= c(-3627271.328111211769283, 1583173.512170388596132),
                            delta= c(1002.008657743770527, 1002.008657743770527),
                            refsys= c(crs_proj4, crs_proj4)
                            )
nc1_out

As my code did not take effect, both nc1_s1 (input to st_set_dimensions()) and nc1_out (its output) show the same dimensions as below.

stars object with 2 dimensions and 2 attributes
attribute(s), summary of first 1e+05 cells:
    X..DQF          X..Rad
 Min.   :  0.0   Min.   : 157.0
 1st Qu.:  0.0   1st Qu.: 256.0
 Median :  0.0   Median : 329.0
 Mean   : 36.5   Mean   : 408.8
 3rd Qu.:  0.0   3rd Qu.: 396.0
 Max.   :255.0   Max.   :1023.0
dimension(s):
  from   to offset delta refsys point values
x    1 5000     NA    NA     NA    NA   NULL [x]
y    1 3000     NA    NA     NA    NA   NULL [y]

How can I set dimensions to NetCDF using stars package?

[Update] This is my environment (see comment by Spacedman that this is not reproducible).

R version 3.5.2 (2018-12-20) -- "Eggshell Igloo"
Running R in Rstudio Version 1.1.463
Windows10 (64bit)
library {stars 0.3-0} is linked to {GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3}


texts:
Is there a way to create a feature layer (or other queryable layer) from a local JSON file? I have a feature layer that was exported to JSON (not geoJSON) and need to display it completely offline. Is this possible?

Using ArcGIS JavaScript API 3.21


texts: Unloading plugin fails after changing toolbar location

texts:
For a Qgis plugin I want to make a toolbar with additional buttons, which I want to dock on the vertically on left side beneath the add layers buttons
I made the toolbar with in the a :

class Test:
    def __init__(self, iface):
        self.iface = iface

    def initGui(self):
        self.toolbar = self.iface.addToolBar("My tools")
        self.iface.addToolBar(self.toolbar, area=Qt.LeftToolBarArea)

To unload the plugin I used:

    def unload(self):
        del self.toolbar

The toolbar is in the proper location, however when I unload the plugin it does not delete the toolbar. And when I reload it it creates an extra instance. So, I appear to have no way to unload the toolbar. I tried assigning a variable to the second addToolBar function but it returns None. Does anyone know of a solution?


texts: Degrees Minutes Seconds to Decimal degrees for Iceland coordinates - Normal formulas come out wrong

texts:
I want to put a file with Coordinates in Minutes seconds to decimal degrees.

Ive tried to use formulas such as should work in excel both modules and this formula: =-H2-I2/60-J2/3600 and 65,26583333. As well as tried to use python in ArcGis. None of them work and the coordinates end up being close to where they chould be but in the sea.

I can see that the coordinates I have are not wrong because Ive checked there status according to minutes seconds format. One coordinate: 65° 03' 839" 22° 44' 982")

Does anybody know what might be wrong or if Iceland in a north location needs a different formula?


texts: Algorithm for splitting polygons/polylines into multiple sections at 180° E/W

texts:
I have code (in C++) that draws various polygons on a map that uses a cylindrical projection (typically Mercator, sometimes Miller Equidistant if you want to be precise).  For simple polygons, curves are drawn between vertices through a set of tiepoints such that the line roughly forms either the great circle segment or rhumb line between the two vertices.  The more complex polygons also include arcs which generally don't form part of a great circle or rhumb line.  At the time of drawing, the tiepoints that make up the polygons are translated into screen positions and drawn using built-in Windows C++ drawing commands.

For the most part, all my polygons and polylines are drawing fine, but I'm encountering problems when at least one line crosses the 180° meridian.  My display is effectively limited to 180° W to 180° E.

One solution I've considered is to unwrap the points that make up the polygon, perform some mathematical wizardry to calculate the screen position for points outside of [-180, 180], and then draw the object twice, once for 180° W and then again for 180° E.  Unfortunately this method falls apart if there are 2 consecutive crossings in the same direction.  As unlikely as that is to happen in my case, the mere possibility has me worried.  It also doesn't work for objects that enclose one of the poles.

Is there an algorithm for splitting a polygon along the 180° meridian while also properly closing it?  Bonus points if it also works for polygons that enclose one of the poles.  Extra bonus points if it'll also work for polar and/or conic projections (not currently supported, but may be required in a later project).


texts: Point values obtained from rasters with different extent, projection and resolution

texts:
I am trying to predict a species' distribution. I am fairly new to GIS, so apologies if I am talking nonsense.
I gathered various geo-referenced environmental data from different databases. The files come in slightly varying extents (global scale), most have WGS84 projection but some do not, and generally the resolution varies between 250m to 1km. I have around 500k occurrence points and intended to sample 500k points randomly as background data. I read that reprojecting rasters should be avoided at all costs and instead the spatial points should be reprojected to the CRS of the particular raster file before extracting values.
I just started reading the introduction to species distribution models by Robert J. Hijmans and Jane Elith who note "For any particular study the layers should all have the same spatial extent, resolution, origin, and projection.". I was a bit surprised to read this as I did not see any harm in extracting point values from different resolutions, extents and projections (as long as points are reprojected).
I hope there are more experienced GIS modelers here that could enlighten me whether I was too naive in simply extracting point values from all these different databases without worrying about harmonizing all the rasters to the same extent, resolution and projection.


texts: Annotations with images in the background

texts:
I was wondering if it is possible to insert text into a QGIS map not just in composer, similar to the Insert Text in ArcMap? I'm trying to put in custom labels with an image as the text box. I've tried the text annotation function but I can't make that work. Is there any other options?

I need to place a text box on a map and it has to be visible there and in composer so I don't have to add the text box to every new composer created.


texts: Creating feature from part of state boundary in shapefile of U.S. states?

texts:
I have a shapefile of U.S. states.  I want to create a feature layer from it so that I can do a Near analysis with it.

Suppose that I am interested in the Georgia boundary, but only the part of Georgia that touches South Carolina?  How could I extract just that segment?

Or, suppose I am additionally interested in the part of Georgia that touches both North Carolina and South Carolina but not at all Alabama, Tennessee, or Florida?


texts: gdal.Open only returns dataset from geotiff when aux.xml file available

texts:
With some specific files, gdal.Open() in python only returns a valid dataset if there's an *.aux.xml file available. The aux file has been generated by QGIS when opening there manually. If no aux file, gdal.Open returns a NoneType.

It is all described in this github repo:
https://github.com/awesomemaptools/gdalissue

Here's the python script:

import sys
from osgeo import gdal
filename = sys.argv[1]
ds = gdal.Open(filename)
print(type(ds))
band = ds.GetRasterBand(1)
data = band.ReadAsArray()
print('value at 1,1:', data[1][1])

When no aux file is available in specific cases:

$ ./gdalOpen.py tmpfile1.tif
ERROR 1: _TIFFVSetField:tmpfile1.tif: Null count for "GeoDoubleParams" (type 12, writecount -1, passcount 1)
<type 'NoneType'>
Traceback (most recent call last):
  File "./gdalOpen.py", line 16, in <module>
    band = ds.GetRasterBand(1)
AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

Any suggestions?


texts:
I have digital orthophotographs from 2013 that I'm attempting to classify in an effort to extract all the cultivated/pasture land from.  I'm using the "Reclassify by Ascii File" tool for this procedure.  According to the NLC database values 81-84 are what I'm looking for and this seems to be correct.  The issue I'm running into is that the classification seems to be the opposite of what I need.

My remap file called "remapfile.txt" is as follows:

81 84 : 1

I understand this as extracting pixels with values 81,82,83, and 84 and remapping them to a value of 1.  I check the box to send all other values to NODATA.  However, what is returned seems to be the opposite.  All the areas that have values of 81-84 are mapped as NODATA and all other points are mapped to a value of 1.

Is there a better way to do what I'm doing???

I'm fairly proficient in Python if that helps.


texts:
I'm new in QGIS and I have been having an issue with it.
The thing is that I'm trying to get the slope and relief of a zone. I have the contour lines, from that I created the DEM with 'Rasterize' > 'Fill no data' and it looks pretty good to me. But when I apply 'Slope' and 'Relief' it creates a weird interaction:
(Slope between 0 and 81,9)

(Zoomed in)

(Relief)

Someone knows what is going on?


texts: When I export a layer with a definition query, it does not maintain coverage. Why is that and how can I find a way to make it do that?

texts:
I am trying to use the NHD waterbody layer in my map. As is, it is a really dense dataset, so I am narrowing it down to only show lakes (Ftype = 390). When I use the definition query, everything looks fine, both in the map viewer and in the layout view. When I export the layout, it appears that suddenly the waterbody coverage expands. It doesn't expand to the full NHD coverage, so clearly some filtering is being maintained through the export, but I'm not sure what is causing this change?
I tried using Select by Attribute (Ftype = 390) and that result seems to match the expanded coverage that results from exporting. To be honest, this is still too dense and distracts from the other data. I like the way it looks when I use the definition query. So two questions:Wwhy are the definition query and select by attribute tools spitting out two different results when the expression is the same? Also, how can I save the result from the definition query?
I am using ArcGIS Pro 2.6


texts: Standard(s) to store and share geolocation and 3D orientation data for an object, e.g. a raster

texts:
As an example, GeoTIFF is a great standard to embed geolocation data into a raster file for cartographic purposes.
But in that case, raster files are usually supposed to be on the projection plane.
What if I want now to share the position (i.e. geolocation) + 3D orientation of an oblique raster (perspective projection) e.g. a simple photograph for example, for which the position and the orientation data are known? For example its position (X,Y,Z) is given as WGS84 coordinates + ellipsoidal height, and its orientation is given as 3 locally defined attitude angles, such as yaw(Z), pitch(X) and roll(Y).
Does there exist any well established and used (geo-)standardized way to store and share these data, either both together position + orientation (best case) or separately?


texts: Point in polygon algorithm for multiple polygons

texts:
I've a Google map with bunch of polygons on it.

Here is a problem I'm interested in: Given a lat,lng point, what's the best way to determine all the polygons that this point lies in?

The obvious way is to run a "point in polygon" algorithm iteratively for each polygon, but I was wondering if there an efficient algorithm to answer such queries especially if you have thousands of polygons.


texts: How do I remove roads classified as residential from an OSM slippy map?

texts:
I'm developing an application in which I need to have an underlying graph representation of streets. Google Maps does not provide this data, which is why I decided to use OpenStreetMap.

I have successfully gotten these networks, but to avoid making my algorithm handle too many edges at a time, I've excluded those marked as "residential" by OpenStreetMap.

Is there any corresponding way I can exclude those same roads in the OSM embed?

I've tried Leaflet and OpenLayers and I cannot find any documentation regarding limiting certain types of information.


texts: I'm attempting to convert a .dwg to .shp with FME software. So far I have pulled in the .dwg, added the writer to convert to a shapefile with a specified coord system and chose the location of the output dataset, and now my screen looks like this but I do not know how to actually access or move/save these shapefiles so I can open in arcmap.


texts:
I'm a quite new user of GeoServer.

I've a created a layergroup, and it's work with the layer preview (OpenLayers).

But if I save the source as a .html to the directory of data dir/www, and open it I can see the layers, but can't get the feature data.

What do I wrong?


texts: Can I calculate an average height from a contour dataset?

texts:
Basically I have a contour dataset and need to find the average height.. seems simple but I cannot seem to find an answer.

Secondly... I am trying to find the average height of selected area within the dataset. so far I have selected the area with the  select tool and viewing selected attributes on the attributes table. Can I calculate the average of the selected contour lines in the attributes table? E.g. calculate average of a selected field?


texts: Combining line vector shapefiles with different attributes

texts:
I'm using ArcGIS 10.3 for desktop.

I have two shapefiles, both with polyline vector data for roads in the same area. This area is quite large, approximately 7,500 km², so there are a lot of roads in both files. One of the datasets contains condition rating information for the roads, the other has the width of the roads. It is necessary for me to combine these two files together so as to have the attributes in the same shapefile to work on.

The width dataset contains more roads and is better in terms of connectivity than the condition dataset. I have tried creating a buffer around the condition data then using the Identity tool to combine that data with the width data that overlays it. However this creates a lot of slivers of road at junctions that need to be corrected manually. While this is okay for a small area it is unworkable if I scale up to the whole county area. Is there a way of fixing this? Or is there a better way to combine the line datasets?


texts: Openlayers panel control editing after initialization

texts:
How do I reference controls that have already been added to a panel.  I am attempting to destroy or modify the default controls associated with the editingToolbar.  These changes would be done outside of the initialization function.

EDIT:My intent was to change the layer that a drawFeature controls was adding features too.


texts:
This is my project:

I want to convert all points that touch each other into polygons, something like this (I know I missed some of them):

I tried using dissolve and convex hull to do that, but they didn't work as intended. Also, this is just a fraction of the whole project, so I don't know how many clusters there should be. Hence, I can't run k-means.
Unfortunately, I don't know anything about PostGIS. Does anyone know how to do this?


texts: Classified image using RF method in R displays less no. of classes than it was trained to

texts:
I am trying to classify crops using Sentinel-2 image including 10 classes, using randomForest classification in R. For that, I extracted 2000 random samples points from the training areas for 11 bands in sentinel-2 images and stored them in a dataframe to train and validate the model. It consists of 11 predictors (bands) and one column containing response variable (Class). The structure of the dataframe looks like this:
Classes ‘data.table’ and 'data.frame':  2000 obs. of  12 variables:
$ Clss: Factor w/ 10 levels "Alfalfa","Barren",..: 4 4 2 5 2 4 5 3 4 1 ...
$ B2  : num  0.122 0.1621 0.2427 0.0267 0.3304 ...
$ B3  : num  0.1812 0.233 0.3134 0.0548 0.4002 ...
$ B4  : num  0.2907 0.356 0.3995 0.0203 0.4682 .
$ B5  : num  0.3319 0.3861 0.4245 0.0672 0.4894 ...
$ B6  : num  0.35 0.4 0.429 0.379 0.48 ...
$ B7  : num  0.382 0.429 0.441 0.554 0.491 ...`
$ B8  : num  0.424 0.495 0.462 0.579 0.507 ...
$ B8A : num  0.435 0.489 0.444 0.577 0.485 ..
$ B9  : num  0.447 0.492 0.452 0.573 0.473 ...
$ B11 : num  0.482 0.516 0.544 0.19 0.574 ...
$ B12 : num  0.3169 0.361 0.5553 0.0844 0.5559 ...
- attr(*, ".internal.selfref")=<externalptr>

Then, I performed 10 fold cross validation and grid search to optimize the model using 'trainControl' function. Following is the code to generate final fitted model.
rfmodel<- train(factor(Clss) ~ . , method ="rf", type="Classification",data=B.train, importance=TRUE,
allowParallel=TRUE,tune.grid=data.frame(mtry= c(2:11)) , ntree=800,
trControl=ctrl,metric="AUC")

After fitting the model, I predicted the whole raster image using the 'predict' function in raster package using the code below:
img_predict <- predict(img,rfmodel,type='raw',progress='window')

Here, img is my raster image, that I am trying to predict and rfmodel is my fitted model.
The output of the predicted raster 'img_predict' looks like this:
class      : RasterLayer
dimensions : 8359, 6583, 55027297  (nrow, ncol, ncell)
resolution : 10, 10  (x, y)
extent     : 602990, 668820, 3608270, 3691860  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs
names      : layer
values     : 2, 9  (min, max)
attributes :
   ID   value
from:  1 Alfalfa
to : 10   Wheat

Now the problem is that I trained the classifier using 10 different classes, however the classified raster has only two classes. Though the attributes in the predicted raster shows 10 classes, the plot of the raster showed only two classes. This seems quite strange to me. Am I missing something while using 'predict' function?
Some of the classes in my training sample contains relatively larger sample points than other. Although, I am getting quite good accuracy for the classifier results, could that be the problem? What should I add or check to get proper classified image with all 10 class?
I've tried this multiple times but getting the similar results and could not identify the reason. I am quite new to R. Any help would be appreciated.
This is the first time I am posting something here, so I apologize for the formatting issues.
RGB composite of the image is shown on left and the predicted raster on the right. As seen below, my input raster is not exactly rectangular. However, R is predicting that outside portion as one class and whole study area as one.


texts: Is it possible to set a tool as optional within ArcGIS Model Builder?

texts:
I have a model with two output files - one is useful always and the other only under certain circumstances. The first (always wanted) is a raster dataset, the other (optional) is a vector feature layer derived from the raster.

Is it possible to make the last step (including the corresponding tool) optional, or do I have to make two separate models, one with and one without the final step?

EDIT: I use ArcGIS 10.1.

The output feature layer can help, but it's hard to interpret, so I usually hid it just after creation. It was originally intended as a first step of subsequent model and it will be moved to it, but there are some problems in the next model's creation. Having this layer for each set of parameters is quite handy for me to work on the subsequent model, but almost useless for other users.

Plus, there's curiosity behind that - quite a lot of tools have optional outputs, so even if I didn't need it I'd be curious whether it's possible in model builder too.


texts: Understanding Cut Polygon tool of ArcGIS for Desktop?

texts:
I try to cut this polygon layer into parts

and didn't succeed to cut this area:

I get this error: "A polygon cut operation could not classify all parts of the polygon as left or right of the cutting line."
I red

Cutting Polygon in Two in ArcMap?
How to cut polygon or split polygon using polyline and ArcObjects?

but didn't find an answer to my problem

But i don't understand why i did succeed to cut this area in the polygon:


texts:
I installed ArcGIS 10.2.2 on my laptop and try to localize the ArcMap component to my native language. I used programs like (Lingobit localizer and SDL passolo) to translate resource strings in the file ArcMap.exe but no effect occurred on ArcMap interface. I noticed that some DLLs related to UI should be translated too, but there is a huge number of them.

Anybody can direct me to know which Files (DLLs) should I translate to make ArcMap user interface in my language?


texts: Strange problems with project crashes and bad:alloc in QGIS 1.8 linux

texts:
during last days I had to prepare some maps from field data. I mainly worked with point shapefiles and openlayers map background.

After some time I always noticed some problem like just partial refresh of the background map (although the connection was fine) or much slower behavior of QGIS.

I do not think the CPU and RAM load was big - there was only one shapefile active and small area from openlayers (lets say about 1x1 km) loaded. About 4 maps in map creator.

I always tried to save the project (looked to be successful) but after QGIS restart it always crashed during loading it or showed bad:alloc message and froze.

First though about problem on my laptop but later moved my data to desktop PC (both Ubuntu 13, KDE4, QGIS 1.8 from QGIS master) and although created a new project from scratch, still had problems mentioned above...

Is there something wrong with the openlayers plugin?

anyone else with similar problem?

thanks


texts:
I have an FME workspace, part of which looks like this:

Incoming features may have point geometry or may not, this is because source data is smallworld, where geometry is a non-mandatory attribute.

The CoordinateExtractor transformer is set up like this:

When incoming feature has no geometry, return value for _x and _y attribute is #1.QNAN which means 'quiet not a number', as I've learned so far, and means something like NULL value.

The question is, how to test for this value in the Tester transformer, as far as I see, Attribute has null value test does not work.


texts: Error arcpy.da.InsertCursor from list to table

texts:
I'd like to export a FieldList into a new table, but I get an error:
Traceback (most recent call last):
File "D:\TOOLBOX\kmap1_mod.py", line 194, in
cursor.insertRow(row)
TypeError: sequence size must match size of the row
Here is my code:
field_names = [f.name for f in arcpy.ListFields(map2, "*T3*")]

tabl=arcpy.CreateTable_management(topo, "datum")
arcpy.AddField_management(tabl, "F1", "TEXT", 60)
with arcpy.da.InsertCursor(tabl, ["F1"]) as cursor:
    for row in field_names:
        cursor.insertRow(row)
    del cursor

I would like to get this result in the final Excel table:

How can I solve this problem?


texts: IElement added to IGraphicsContainer ignores zOrder

texts:
I display a balloon callout after casting a axMapControl.Map to a IActiveView and doing this:

pGraphicsContainer = activeView as IGraphicsContainer;
pGraphicsContainer.AddElement(pElement, zOrder);
pElement.Activate(activeView.ScreenDisplay);

The problem is the balloon callout displays below the place labels layer. The text "Minneapolis" is visible on top of my balloon callout.  I have tested zOrder with values from -10 to 10000 and the place labels still render on top of the balloon callout.


texts: Exporting separate files of streets from MID/MID which cover district

texts:
I have MID/MIF that describes 80 streets of the district and I want to have separate pairs of MID/MID for each street from the original MID/MID.
I can save a copy of district MID/MID and delete all but one street and save it, but maybe it could be made automatically?
Here is original MID/MIF:

so I need separate MID/MIF for each table entry.


texts: Seeking up-to-date shapefile(s) for Pakistan?

texts:
I'm currently looking for a shapefile for Pakistan with the current district borders. So far I found something on https://gadm.org/. Unfortunately the shapefile there seems not to be up-to-date, as e.g. one district has been split into two in 1996, but in the shapefile it's still one district (and in some other cases the names seem to be the old ones).

Is there any "general source" with newer shapefiles (also for other countries)?


texts: Displaying all townships a lake falls in

texts:
I have a shapefile for lakes, and one for townships. Both are polygon files and all lakes fall inside at least one township and I'm using ArcMap. I'm trying to find a way to quickly and simply generate a table that would list the lakes and then all the townships that lake falls in.
Ex:

Lake       Townships
Lake1      Township1, Township2, Township3

I unioned the lakes with the townships to get a new file that had the lakes and townships but the records for each lake were broken up, I need a way to get all the townships into that one single field.


texts: ArcSDE error: "Maximum number of connections to instance exceeded "

texts:
I have an ArcSDE error: " **Failed to connect to database. Maximum number of connections to instance exceeded "*

Did what ESRI recommended :"

  Check the value set for the CONNECTIONS parameter in the SERVER_CONFIG
  or sde_server_config ArcSDE system table to be sure it is sufficient.
  If it is not, increase the number of allowed connections using the
  sdeconfig command with the alter operation"

and still nothing works on my database server. (I can't create features, import, etc).

We use a Microsoft SQL Server. And it's an enterprise geodatabase.

Apparently it was a schema lock. I had a map service started which used that geodatabase. So how can I edit the geodatabase and have a started map service in the same time. I've already created several versions and tried to publish just the versions, but they still lock my gdb.


texts: Rounding raster values in QGIS raster calculator?

texts:
Is it possible to round raster value in QGIS to a given number of digits? I have raster data with temperatures as values but the authors have made the values Celsius * 10. That is, a temperature of 2.5 is stored as 25. When I use the following formula in the Raster Calculator:

"july_merged@1" * .1

...the result is 2.599999904632568. I just want the 2.5.


texts:
I am running python in the ArcGIS python window.
ArcGIS shows error message when trying data merge from arcpy.

File "", line 1, in    File "c:\program files
(x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 3816, in
Merge     raise e ExecuteError: ERROR 000732:

import os
import arcpy
arcpy.env.workspace = "C:/the shapefile folder/"

path = os.listdir("C:/the shapefile folder/")#unzipped all shapefiles to this folder
shapefile = []
for i in path:
    suffix = ".shp";
    if i.endswith(suffix):#getting all files end with '.shp' to a new list
        shapefile.append(i)
arcpy.Merge_management(shapefile,"C:/output folder")

Where is the problem?


texts: Create map of frost free days using Daymet data in ArcGIS

texts:
I stumbled on DAYMET today for creating graphs in R for frost free days at a given location. This was all using boiler plate code someone in my agency made where you just change lat/long; I'm an R and python idiot.

I was wondering if there is a way within Arcmap itself where I can get the consecutive frost free days for each year out of the netCDF files, then create contour lines or just a shapefile with the different bands of FFDs averaged out from 1980-2018. When I search similar ideas, most of them involve R, Python, or google earth engine. Is there an Arc exclusive way to do this?


texts:
I've been running Excel to Google Earth Dashboards for some time through VBA on Windows 7 and Excel 2010. However, I have just tried to copy one of my workbooks to a Windows 8.1 machine on Excel 2016 and the Google Earth 1.0 Type Library is missing from VBA reference so the Google Earth links will not load.The same version of Google Earth is being run on both machines.

I have also tried this on a Windows 10 machine and cannot get it to work either. Any ideas why the Google Earth 1.0 Type Library is not available and where I can get it? A search online returns little information.


texts: Retrieve attribute with max value by group

texts:
I have a intersected counties shapefile with wind potential category attribute. What I want to do is to dissolve the features grouping them by county, and retrieve the WPC value which has the largest AREA. With Dissolve tool i am able to retrieve the area only.

Here is the data snippet:

So my desired output would be a layer with 1 feature for county with WPC value (with max Shape area).


texts: how do i rotate a polygon drawn from the Create Features menu?


texts: Force Arcmap to respect primary key constraints/triggers

texts:
We have a feature class that has a primary key constraint(not OBJECTID) on it's database table. This field also has an auto incrementing trigger (BEFORE INSERT)

When our users try to create/add features to this class using arcmap, since it's a non nullable number field, ArcMap adds a default value to it (0), so the second time it tries to put 0, it violates the pk constraint. Not only that but it overwrites the BEFORE INSERT trigger by putting 0 after it triggered.

I haven't been able to workaround this even with indexes.

An idea that I thought could work was making a AFTER INSERT incrementing trigger, or an ON UPDATE trigger that only works if the pk values is 0, but I'm not sure how to make none of these, or if that's even possible/necessary.


texts: Speeding up Python code to convert multiple csv to shapefile?

texts:
I am trying to work with some crime data from https://data.police.uk/data/
The data is organised in several .csv files, one for each month and each crime is geocoded with Lat and Long.

As the file structure might differ from month to month I cannot merge the csv together using
copy *.csv combined.csv in the command prompt as explained here: https://www.itsupportguides.com/office-2010/how-to-merge-multiple-csv-files-into-one-csv-file-using-cmd/

So I decided to use python to loop through all the csv files in the folder and create a shapefile for each one which I will then merge together at a latter stage.

This is the code I came up with after looking at this post Convertion of multiple csv automatically to shp, it works but it is really slow, in a couple of hours it converted only a handful of tables. Do you have any suggestion to improve my code?

I had to use csvfile.replace('-', '_') as the file names look like this 2012-05-metropolitan-street.csv and I cannot use "-" in the output shapefile name.

import arcpy,os
shpworkspace = r"G:\GIS DATA\Crime Data\CSV"
arcpy.env.workspace = shpworkspace
arcpy.env.overwriteOutput = True

csvlist = arcpy.ListFiles("*.csv")

try:
    for csvfile in csvlist:
        outlayer = "CSVEventLayer"
        spatialreference = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision"
        arcpy.MakeXYEventLayer_management(csvfile,"Longitude","Latitude",outlayer,spatialreference,"#")

        shpfile = os.path.splitext(csvfile.replace('-', '_'))[0]
        arcpy.CopyFeatures_management(outlayer,shpfile)
    del outlayer

except:
    # If an error occurred print the message to the screen
    print arcpy.GetMessages()


texts: Generating Scale Bar for Raster Image Data

texts:
My goal is to create a scale bar for a printed map at 1:1,000,000 scale.

The input dataset is a Mercator projected geotiff.  Using GDAL it is possible to get the spatial reference information and dataset metadata.  One component of the metadata is pixel resolution (at the equator).  Assume each pixel is ~4km square.

The scale bar will be simplistic, and the goal is determine the printed length of the scale bar at 100km.

I am struggling to define the ratio between printed size and visualized size.  For example, at 1:1000000, I can define the unit as cm and know that a 10cm scale bar is required.  But this needs to be related back to the map and I am completely missing that (basic) linkage.

I am simplifying and assuming the scale bar only needs to be accurate at the equator.


texts: Join attribute by location result lacks some features in QGIS 2.18

texts:
I'm new in QGIS and I'm trying to do a Join attribute by location and  Count points in polygon using a polygon shapefile and a points shapefile. After running these tools, an error message is printed:

"QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType' See log for more details"

Even tho it creates a shapefile, I look into the returning shapefile and it creates some features over the area where points are located but it doesn't create features where other points are located. It would be great to know the cause of this problem and how to fix it.


texts:
So I have a script that calculates a diameter of a crater of an asteroid with user input parameters. what i need help with is creating a buffer that brings in the diameter as the distance value and creates a buffer around different point locations. I can't seem to find a way to bring the value into the buffer tool. I have also tried to find a way to just add the distance to layer i want to do a buffer around as an attribute.

def main():
    import arcpy
    velocity = arcpy.GetParameter(0)
    mass = arcpy.GetParameter(1)
    component = GetParameterAsText(2)

    dct_mat = {'ice' : 1.8, 'iron' : 7.3, 'rock' : 2.3}
    comp = dct_mat[component]
    grav = (0.16667)
    const = (0.007)
    colrt = (1.3)
    soil = (2.65)
    power = (1.0/3.4)
    vel = (velocity * 1000.0)
    kilotons = (4184000000000.0)
    kinetic = (0.5 * mass * (vel ** 2))
    energy = (kinetic * kilotons)
    part1 = (const * colrt * grav)
    part2 = (comp / soil)
    part3 = (energy * part2)
    part4 = (part3 ** power)
    crater = (part1 * part4)
    final = (crater / int(1000))

    arcpy.SetParameter(3, final)
    arcpy.AddMessage("final: {}".format(final))

if __name__ == '_main_':
    main()

here is the script i'm using. also i cant get the message to work but that isn't as important right now.


texts:
I have a database table with a dynamic schema implemented via a JSON column named "properties".

Currently GeoServer doesn't recognize my JSON column.

Is there some way to make Geoserver recognize the column?

Or may be if it doesn't work with JSON, is there some way to make it work with HSTore?

I need the properties JSON column to be rendered in the properties of a WFS response, also I need to apply several dynamic CQL_FILTER to the "properties" column.


texts: Gdalwarp UTM reprojection causes image distortions

texts:
I'm trying to re-project Sentinel-2 images to a common UTM zone in order to create a mosaic.

For this I'm using gdalwarp with the following syntax:

C:/OSGeo4W64/OSGeo4W.bat; gdalwarp -t_srs \"+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0\" -tr 10 10 img1_UTM35s.tif img2_UTM36s_reproj.tif

While the execution works fine, the results show some funny distortions. It seems like groups of pixels inside the image rotate clockwise or counter-clockwise which leads to a small shift as seen in this image:

The right side shows the same part of the original image as reprojected to UTM36S. Inside the red circles you can see (albeit barely) the shift of the pixels resulting in apparent edges.

Did anyone experience this already?


texts: Insert points to PostGIS using a function

texts:
I am trying to insert points geometry into PostGIS database. Well I tried to use a function just to read point coordinates and attributes and save them in a table in PostGIS. It works for attributes but not for geometry. This is my python code:

insert(id, Height, X, Y)
....
def insert(gid,height, X, Y):
   gid=str(gid)
   connection = psycopg2.connect(host='', dbname='postgres',   user='', password='')
    cur = connection.cursor()
   cur.execute ("""INSERT INTO gronkarta.test5 (gid, height, geom) VALUES (""" + gid + """,""" + str(height) + ""","""+ST_GeomFromText('POINT(X Y)', 3008)+""");""")


texts: How aggregate data in a geodataframe by the geometry in a geoseries

texts:
I have a geoseries where each row is a (square) polygon. Housing price data is in a geodataframe where the geometry is a polygon of the zip code.

I need to aggregate the housing price data by squares in the geoseries.

For example for a square that is within a zip code the value would be the housing price of the zip code. When the square in the geoseries intersects many zip codes I would like to get the weighted average of the housing prices in the zip codes, where the weights would be the shares of the areas intersected by each zip code. (Or if this weighted average is difficult, I could take the price from the zip code whose share is greatest in the square)


texts: how to highlight centroid of a radius area when zooming to extent of radii

texts:
I am trying to figure out a way to highlight the centroid of a radius after right clicking on the radii attribute and zooming to its extent.  ultimately, I want to batch export a few hundred maps with each map showing the highlighted centroid at a 5 mile extent.
Another thought I have (if its possible) is to set the 'zoom to selected feature' extent to 5 miles for the points shpfile.  Then I can discard the radius shpfile and work soley with the points.  Is there a way to do this?


texts: Copying multiple .shp files to a file geodatabase

texts:
I am trying to copy multiple shape files to a file geodatabase using python. I have been able to copy single .shp files to a geodatabase. Doing this one at a time is tedious and seems easier to do in ArcMap. I have enclosed the code I am using to do one at a time. How can I make this script work for all .shp files in the folder to be copied to a single geodatabase within that folder?

# Name: CopyFeatures_Example2.py
# Description: Convert all shapefiles in a folder to geodatabase feature
# classes

# Import system modules
import arcpy
import os

# Set environment settings
arcpy.env.workspace = "C:/gisdata/roadData"

# Set local variables
outWorkspace = "C:/GISDATA/RoadData/output.gdb"

# Use ListFeatureClasses to generate a list of shapefiles in the
#  workspace shown above.
fcList = arcpy.ListFeatureClasses()

# Execute CopyFeatures for each input shapefile
for shapefile in fcList:
    # Determine the new output feature class path and name
    outFeatureClass = os.path.join(outWorkspace, shapefile.strip(".shp"))
    arcpy.CopyFeatures_management(shapefile, outFeatureClass)


texts: Upgrade strategy for geoserver 2.1.4 to 2.5.x

texts:
I have a 2.1.4. geoserver installation and I wish to upgrade to the latest version which is currently 2.5.2. Can I do this at once by installing the latest version or is it wise to do an intermediate step via geoserver 2.2. and  upgrade further from there?
I noticed for instance that since geoserver 2.2.0 some security changes have been applied. Is it wise to get past this upgrade first and then upgrade further?


texts: Trying to get Multiple Tabs to display in Attributes Table in CMV viewer javascript file

texts:
My first tab displays correctly..

This is my code in my viewer:

Both the service URLS display the attributes table in ArcMap correctly and both the URLS work in a browser

This is the error message i get from growler when attempting to load the second tab the search seems to fail..

I get a Search Error msg that displays saying "Sorry, an error occurred and your search could not be completed." then the message disappears....


texts: Understanding ArcGIS Map Package Details?

texts:
Using ArcGIS 10.6.1 to create Map Package, I just notices Esri is creating two folders v10 and v106 which both holds same .mxd and Geodatabases

It is also different when Runtime enabled and I am getting only v101. couldn't find any documentation on Esri for this so what are these v named folders and ...

Can you please let me know why two files are generating?


texts: Creating Human Trafficking Heatmap in QGIS?

texts:
I'm trying to create a heatmap with about 325 points distributed around the city of Houston. (FYI Houston is a major human trafficking in the US: a title shared with LA and Atlanta depending upon seasons and large events).

I'd like to plot the concentrations of these points, but can only get a large circular area around all of them with QGIS' heatmap plugin.

Currently, QGIS only allows whole numbers as its buffer ratio, so I would need to go smaller or possibly correct my CRS. The confirmed human trafficking addresses were geocoded with Google's API, and I have unsuccessfully tried in WGS84, WGS84/Mercator, and Google Mercator.

Is there a better way to do this? Here's a screen cap to illustrate. The results were the same for the aforementioned CRS's.

PS. I'm also already aware that the current heatmap plugin can only create a raster layer with projected coordinates if the buffer ratio is in map units


texts: Add an image to a popup in a Leaflet map

texts:
I'm trying to add an image to an Easy Button activated popup. I understand that anything within () after .setContent is html formatted, though anything I add other than text does not work. Also when I delete everything within the () and just add an image the entire map vanishes.

code:

  var helloPopup = L.popup().setContent("Portland Cycle Safety Map
This map is intended to bring dangerous intersections
and street segments to the attention of Portland area cyclists: this is a work in progress.
Each skull marks the location a cyclist has been killed by an automobile sometime between 2005 and 2017.");

L.easyButton('fa-crosshairs', function(btn, map){
    helloPopup.setLatLng(map.getCenter()).openOn(map);

screenshot:


texts: Getting error Basemap Layers cannot be published directly to service

texts:
I am trying to take a map from a PDF document and get it so I can use it in ArcGIS Online. I was following a great tutorial but ran into a snag.
My overall methodology is to georeference the pdf map (that I converted into a JPG before putting it into ArcMap Desktop) and then I get it the way I want it. After I georeference the JPG image I rectify and delete the base map I was using. Now I go to share as a service ---> publish a service --> choose a connection ---> and then I fill out the service editor. However, once I click analyze I get this error:
Severity: High
Status: Unresolved
Code: 00037
Description: Basemap Layers cannot be published directly to a service
Name: Basemap
Type: Payer
Data Frame: Layers

I don't know why I am getting this error! I uploaded an image so you can see my screen.
Any advice?


texts: How to export a GTiff with its GDAL colour table, in GRASS, which is suitable for Adobe Illustrator?

texts:
I've created a new raster in GRASS, starting from a dem and using the commands r.relief and r.shade.

Now I would like to export this new raster as a georeferenced image and import it in Adobe Illustrator (using MAPublisher plugin).

I'm trying to do it using the command r.out.gdal and exporting it as GTiff with its colour table using UInt16 data type, with no success...

Have you got any suggestion?


texts:
How do I get the pop-up help window that's part of the auto-complete to show more definitions of the tool being used?

For example in the attached image the definitions of "Buffer_analysis" inputs stop at "line_side".

I'd like to see the definitions for the other inputs as well (Lined_end_type, dissolve, etc...).


texts:
I'm new to ModelBuilder and have successfully used iterators for other things, but I'm having an issue with one particular model...

I have one raster file showing direction to nearest coast ("CoastDir") and I need to reclassify this into 8 separate binary maps representing the 8 compass directions. I figured the easiest way to do this (especially as it's something I have to do fairly regularly) was by creating a geodatabase containing 8 remap tables representing the directions, called N, NE, E, SE, S,... etc.

For example: Table "N" is

My model looks like the image below:

My iterator is set up with the basic geodatabase input and no other settings.

My problem is that it's not iterating properly. It's carrying out the reclassification dozens of times, but only on the "N" table, and overwriting the output, instead of moving onto the next table.

I tried renaming the Input Remap Table as "%Name%" instead of "N", but that didn't work either, but the format above has worked successfully for raster iterations.

I'd rather not delve into python if possible, as I'm building a series of tools, all in ModelBuilder, and I'm even less familiar with scripting. I will also need to do similar reclassifications on more directional rasters.

Do you know what I may be doing wrong?

I've googled an awful lot and am none the wiser! I'm using ArcGIS 10.4.1.


texts: Exporting an animation (Time Series) map through QGIS2web plugin

texts: I simply would like to know if it is possible to export an animation map using the QGIS2web plugin so I may display a time series map on the web through the plugin?


texts: Merge FeatureCollections excluding overlap in Google Earth Engine

texts:
I have 2 FeatureCollections (FCs), roi and fc, that I want to merge. Naively I did it this like below but now I am worried about the overlap in the mergedFC.

var mergedFC = roi.merge(fc);
Map.addLayer(mergedFC, {}, 'merged flooded and non');

The problem is that I don't want my resulting FC to contain duplicates where there is spatial overlap (fc is spatially a subset of roi), I want to keep only the data for fc. This is important because later when I create a training set with .sampleRegions, I'm concerned about functions getting confused about class labels (all roi is non-flooded, whereas all of fc is flooded and of course the latter is what should take priority in an overlap situation).

I thought if I map over fc, taking the .symmetricDifference between each Feature and the roi Feature (there is only one simple Geometry in that FC), then I will end up with a FC that excludes fc because it is a spatial subset of roi. See code below. I thought I could then merge this output with fc which would give me what I want. But when I visualise this output, the whole roi/area seems to be shaded, which suggests it hasn't worked as I expected.

var trying = function(feature) {
  return feature.symmetricDifference(roi.first());
}

var output = fc.map(trying);
Map.addLayer(output, {color: 'FF0000'}, 'output', true, 0.1);

Questions:

Am I correct in assuming there is a problem here, that spatial overlap will introduce class label confusion in a mergedFC that is used for building a training set?
How can I test whether using .symmetricDifference worked?
Is there a better way?
What else might I not have considered?

For context, the example I am learning from is from the TF and EE demo script from the Earth Engine Users Summit 2018.

I'm new to both GEE and to EO.


texts: TypeError: coercing to Unicode: need string or buffer, NoneType found in ArcGIS labeling?

texts:
I have found out how to label using multiple attributes via an expression and python. What I have is as follows:

[ID] + '\n' + [Notes]

An ID example is "MH108" and a notes example is "buried, no access". I have used similar code successfully in the past, but in the current situation I get the following error and cannot figure out what I am missing.

The expression contains and error.
Modify the expression and try again.

Error 0 on line 0
Error running expressions:
Esri_FindLabel(ESRIExpressionArg0,ESRIExpressionARG1)
Traceback (most recent call last):
Files “<expression>”, line 1, in <module>
File “<string>”, line 2, in esri_FindLabel
TypeError: coercing to Unicode: need string or buffer, NoneType found


texts:
All,

Back in December it seemed there was no way to bring in a CartoDB layer into a Mapbox.js or Tilemill project.

Has this changed?

(http://support.mapbox.com/discussions/mapbox-tiles/763-cartodb-and-mapboxjs-how-do-they-talk)

I'd like to implement this on a few local government projects I'll be working on shortly.

Thanks!

-mb


texts: Adding new columns from other shapefiles - Land use density

texts:
New here but I've been following the stack exchange for a while to get many answers on my newbie journey in GIS. I'm trying to find out how to merge different point in polygon classes by column.

Anyhow, I have a land use map of a neighborhood in Chicago. I used the parcel information to create a point map showing all the different land use classes (0-4). From there I went on to create 5 different layers, each showing the parcels of a certain type using new layer by selected features.

Next I did a points in polygon spatial analysis to let me know how how much each land use type exists in a certain block. So I have 5 layers showing me the point count per block - each layer shows a certain type. What I'd like to do is merge the columns. I know I can probably do it by editing the csv/dbf files associated with the shapefile, but I'd like to know if there's a way to do this through QGIS. Also, if you feel like I was inefficient in the way I went about this, please don't hesitate to critique my methodology.


texts:
Programmers' list has been a pleasure to read. Stats' list is growing fast as well. So perhaps there is a place for separate one on this site?

Please share your nuggets of wisdom / humour / inspiration from the world of GIS, cartography and mapping.


texts: How do I create a clickable map where clicking on a specific part of the map takes you to a webpage using QGIS?

texts:
I have a shapefile of regions in Italy. I am creating a travel blog and I want to create a clickable map where when you click one of the regions it takes you to the relevant page for that region. Is there anyway to do this using QGIS? I am also using wix to create the website and I can code in Python if that helps.
Something similar to what is shown here is what I need. You can see for example if you click on the US it takes you to another page.


texts: Seeking free/open source Western Australian Topographical data?

texts:
Does anyone know (or have m) a free or open source Western Australian Topographical datasets?

I have a project that involves terrain relief in Northam, ideally they want it at a scale of 1;25,000 but i thought i would ask around if anyone knew (or had) a dataset before buying one from landgate.


texts: QGIS- Fixing adjacent lines that looks like touching each other but do not

texts:
I'm trying to dissolve and then multipart to single part a line layer. The thing is that after I'm breaking it to single part some lines that look like touching each other break back to little parts. I was thinking maybe they are not really touching each other.
How can I check if they are really touching and, if they are not, how can I fix it?


texts: Create Thiessen polygons using spatially weighted averages?

texts:
I have to create Thiessen polygons with the spatially weighted averages for specific contaminate concentrations from different borings. I have the Thiessen polygons created but this only shows me the polygons based on the point locations not actually point concentration values that I would like. What tool can I use to get this info?

I only have ArcGIS free trial which limits the amount of tools I have access too. I also have QGIS and could use either to do solve the problem.


texts:
I had previously used ENVI to create ROIs for a Random Forest classification. After creating my files in ENVI, I exported these to EnMap and ran the classification there. (Please note knowledge of EnMap is not required to answer this question)

To rerun my results, I created tifs of the satellite data (SPOT 6) & training & test datasets, as this is required for EnMap - or so I thought. I tried to use EnMap again to run my results however I keep getting error messages. After much detective work, I have concluded the reason the software was reluctant to run the classification was because files are not in in ENVI format and tif and bsq just will not do.

I no longer have access to ENVI so cannot convert my files to the ENVI file type for use in EnMap. I do not currently know anyone else with the software who could help. I'm therefore looking for an alternative way to run an RF classification.

I have access to ArcMap and QGIS and am fairly competent with both. I'd be willing to download new software to try however I live in a remote area with limited internet and many freeware cannot be easily learnt in a short space of time, especially when there is no one around to help teach it, as is the case for my situation. So I'm not sure R is an option. I'd still be ok with using EnMap if I could get it to work.


texts: Clip GAPAnalysis and save as a new file output strange looking colors

texts:
I download USGS GAP Analysis and use a polygon to clip raster from that GAP Analysis. I am wondering why there is a strange looking color when you clipped the raster. If you zoom in, you can see myraid colors, but when zoom out you see some white spotty areas.

I can't insert 2 saved image for you to review it, but here is the link you can see it .

Gap Analysis example


texts:
I have a line that goes through sections of land in which has attribute data assigned to it.  I want to compute out how long of pipe goes through each different type of land.
As you can see, the cyan lines are the boundaries of land and this one happens to have the attribute "Cultivated_Crops".

The best way I know of doing this is through autocad just measuring the distance by hand.  I'm not sure of any quick solution.


texts: Why can't I clip vectors with 2 different projections?

texts:
I am using QGIS.  I am trying to clip parcel data with a stream buffer.  The two have different projections.  Parcels are NAD83 State plane feet and streams are NAD83 State plane meters.  I keep getting an error that says something like

detected different projections,  some data may not disply correctly.

When it runs the clip the output file has no data in it. I have enabled on the fly projection.


texts: Where to download ASTER scene boundaries?

texts: I am looking for a shapefile of ASTER scene boundaries for the Kansas and Oklahoma region.  The attached screenshot from GLOVIS is what I'm after, except I need polygon footprints with relevant attributes such as scene ID and year.  Is such a product available for download?


texts: Rasterize Function in QGIS 2.18.13 Gives an Error

texts:
I just updated QGIS to 2.18.13 and now I'm having trouble using the GDAL [Conversion] function, Rasterize (vector to raster).

I have sixteen vectors I made and I want to convert them to rasters. I've done this before with little to no issue, and I have never gotten the error I keep getting.

The message log says:

  'NoneType' object has no attribute '__getitem__'

but it says this before I can even open a window. Literally, I click "Rasterize (vector to raster)" and this error comes up before I can open a new window to put in the vector I want to rasterize, or the attribute that I want to fill the raster. Is this a bug or am I missing something?