texts:
I created a QFormLayout and I want to save the text that user types from the QLineEdit into some variables. So in my code, I want to get user input from e1,e2, e3, and e4. Here is my code:
import sys

class Dialog(QDialog):
    def __init__(self):
        super(Dialog, self).__init__()
        self.createFormGroupBox()

        buttonBox = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel)
        buttonBox.accepted.connect(self.accept)
        buttonBox.rejected.connect(self.reject)

        mainLayout = QVBoxLayout()
        mainLayout.addWidget(self.formGroupBox)
        mainLayout.addWidget(buttonBox)
        self.setLayout(mainLayout)
        self.setWindowTitle("Form Layout - pythonspot.com")

    def createFormGroupBox(self):
        self.formGroupBox = QGroupBox("Form layout")
        self.layout = QFormLayout()
        self.e1 = QLineEdit()
        self.e2 = QLineEdit()
        self.e3 = QLineEdit()
        self.e4 = QLineEdit()
        self.layout.addRow(QLabel("Date debut:"), self.e1)
        self.layout.addRow(QLabel("Heure debut:"), self.e2)
        self.layout.addRow(QLabel("Date fin:"), self.e3)
        self.layout.addRow(QLabel("heure fin:"), self.e4)
        self.formGroupBox.setLayout(self.layout)

app = QApplication(sys.argv)
dialog = Dialog()
dialog.show()
#sys.exit(dialog.exec_())


texts: Way to speed up TIN interpolation (Delaunay triangulation) QGIS 3

texts:
Are there ways to improve processing time for TIN interpolation?

I'm using TIN interpolation as a step in a larger process.  The interpolation step is taking a very long time (10+ hours to a few days), and I often end up with missing areas (see example below).
I'm only using ~2000 input points.  My output raster is a 1m elevation raster grid covering about 25km x 5km so ~125 million cells.  Output raster needs to match the resolution of another layer, so I can't just reduce the resolution.

The input values are essentially coming from elevation contours along a valley slope, but the valleys I'm working with are fairly steep, with frequent gradient changes - there is some room to use sparser input data, but I suspect not enough to cut processing time down to a few hours.

I've tried dividing the data into smaller blocks, but this doesn't seem to save any real time over processing it as a block. It does help with the missing data issues, however.

I've tried saving the output as both memory layers and on disk as a TIF, this doesn't seem to affect the processing time.

I'm using QGIS 3.8.3, have also tried on the LTR (3.4.12).  I'm also open to a solution in R, as I use it for other steps in the process, but haven't found a TIN interpolation algorithm in R.  I prefer the TIN algorithm to IDW or kriging because of how it respects the input contours and goes through known points.

Image with output raster missing an area. I assume that the missing area is the result of my computer running out of resources during the processing.


texts: Improving performance of looping over MXDs using ArcPy?

texts:
I am attempting to find the fastest method for looping over multiple map documents while keeping memory resources down. I'm running this script outside of Arcmap 10.

I have 5 map documents which are about 1mb each.

I have used both return functions and generator functions. After some testing I found the following code below to be the fastest, running in about 20 seconds.

Memory usage balloons to about 220 000 k when running then back down to 150 000 k after the del statements. If I remove the the del statements run in 18 secs but memory after remains at 220 000 k. which is undesirable.

Note: Importing arcpy on its own uses 117 000 k

I am looking for any feedback you have or suggestions to get increase speed, and reduce memory usage for the code below.

import time
import glob
import arcpy

start = time.time()

path = "F:\\Jobs2012\\06292\\ESRI\\"

def getMxdPath(path):
    """ Returns list of mxd filepath
   as strings for a given path."""
    mxdPaths = list()
    for mxd in glob.glob(path + r"*.mxd"):
        mxdPaths.append(mxd)
    return mxdPaths

def getMxdObjects(mxdPaths):
    """ Returns list of mxd objects"""
    mxdObjects = list()
    for mxd in mxdPaths:
        mxdObjects.append(arcpy.mapping.MapDocument(mxd))
    return mxdObjects

def genr(mxdObjects):
    for mxd in mxdObjects:
        yield mxd

def saveMxd(mxdGen):
    for mxd in mxdGen:
        print mxd.filePath
        mxd.save()
        #del mxd

mxdPaths = getMxdPath(path)

mxdObject = getMxdObjects(mxdPaths)
#del mxdPaths

mxdGen = genr(mxdObject)
#del mxdObject

saveMxd(mxdGen)
#del mxdGen

del mxdPaths
del mxdObject
del mxdGen

elapsed = (time.time() - start)

print elapsed


texts:
I have a QGis project that I programmatically modify using the Python xml.etree library.
Now, I would like to know if it is possible to add few lines to my script in order to export the modified project to a new jpg image. The python script should be run OUTSIDE the QGis python console. Thanks!


texts: ArcGIS JS API Legend Dijit not refreshing correctly

texts:
I have a web application built with v3.5 of the JS API. I have some Jquery functionality that changes the map being built based on user input. My problem is that I can't get the legend to stay in sync with the map. On the first selection I build a map with new esri.Map and I create a new esri.dijit.Legend on the map's onLayersAddResult event. With the first selection everything works and is pretty snappy. When a user chooses the next map setup I've tried two different techniques:

The first is to build a brand new map. I call map.destroy() on the old one to get rid of it. I create a new map and add my layers to it but a new legend does not get created on the onLayersAddResult. The event fires but it seems it cannot create a new legend. If i use legend.destroy() before creating a new one, it blows away the whole div and can't create the legend either.
The second technique is not to create a new map. I remove all the layers from the map and add the new layers to the map. Then I point legend.map to the map and use legend.refresh() with the new layers as an argument. This works, sort of. It seems legend.refresh is really slow and if my second map has more than 3 or 4 layers in it (which it does) it will never draw. In fact it hijacks the whole app and everything is frozen.

So can anyone give me some tips on how to use the legend dijit correctly? Either get rid of it properly so I can build a new one or have it refresh correctly and promptly with new map layers?


texts: remove one image from imageCollection in google earth engine

texts:
I work with google earth engine on an imageCollection from landsat7. Some of those images are 'corrupted'. So I would like to remove one image from an image collection.

Actually I have :

var collectionLS7NDVI0 = ee.ImageCollection('LANDSAT/LE7_L1T_32DAY_NDVI')
    .filter(ee.Filter.calendarRange(6,6,'month'))
    .map(addTime);

var listOfImages = collectionLS7NDVI0.toList(collectionLS7NDVI0.size());
print(listOfImages)
0:
Image LANDSAT/LE7_L1T_32DAY_NDVI/19990610 (2 bands)
1:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20000609 (2 bands)
2:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20010610 (2 bands)
3:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20020610 (2 bands)
4:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20030610 (2 bands)
5:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20040609 (2 bands)
6:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20050610 (2 bands)
7:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20060610 (2 bands)
8:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20070610 (2 bands)
9:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20080609 (2 bands)
10:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20090610 (2 bands)
11:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20100610 (2 bands)
12:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20110610 (2 bands)
13:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20120609 (2 bands)
14:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20130610 (2 bands)
15:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20140610 (2 bands)
16:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20150610 (2 bands)
17:
Image LANDSAT/LE7_L1T_32DAY_NDVI/20160609 (2 bands)

And I want to remove images 6 (2005).

Any ideas are welcome

Etienne


texts:
I try to host an OSM map locally on my server. I follow this tutorial Rendering OSM on ubuntu and managed to installed MapServer and MapCache.

The problem is, when I try to zoom the map it can only zoom once and other will just repeatedly create the same tile, sorry I couldn't post the image, since I'm new here.

What I understand is when I enter: http://localhost/mapcache/demo/tms the Apache web server will look for XML file (mapcache-osm.xml) and the XML file will find call the Mapfile and then it will generate the tile.

Is it true? Please correct me.

I still couldn't figure out how to make the zoom right, did I miss something on installation?

I also test the TMS on OpenLayers site as recommended by the tutorial site. But the result is the same.

I thought if the Mapfile and database (PostGIS) connect successfully then it will generate the tile and I will be able to see all the tile at any zoom level. Or is there something that I need more to create the zoom features?

I found that if I enter http://localhost/mapcache/demo/tms then the XML file generate openlayers automatically. Am I wrong? Please correct me with the concept.


texts: GeoServer features in layers are not displayed in GML when requested by their FeatureID

texts:
My GeoServer WFS service behaves strangely:

http://myserverurl.de/geoserver/my:layer/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Ahaus_neu:Graeber&maxFeatures=50&

I get 50 features, as I expect, and the features are tagged like this:

<gml:featureMember>
  <my:layer fid="myfeature.1">
    ...
  </my:layer>
 </gml:featureMember>
....
<gml:featureMember>
  <my:layer fid="myfeature.50">
    ...
  </my:layer>
 </gml:featureMember>

Usually, I could retrieve single Feature myfeature.XY like this

http://myserverurl.de/geoserver/my:layer/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Ahaus_neu:Graeber&featureID=myfeature.XY&

This works with featureID=myfeature.1, featureID=myfeature.2, and, weirdly, featureID=myfeature.387

It does NOT work with any other featureID, even though they are displayed in the main request.

For featureID=myfeature.3, for example, I just get:

<wfs:FeatureCollection xsi:schemaLocation="myurl//WFS-basic.xsd">
   <gml:boundedBy>
     <gml:null>unknown</gml:null>
   </gml:boundedBy>
</wfs:FeatureCollection>

There is nothing that featureID myfeature.1, myfeature.2, and myfeature.387 have in common except that. This behaviour correlates with the ability of any property of the feature to be updated from OpenLayers. (Which works only for myfeature.1, myfeature.2 and myfeature.387).

Any ideas?


texts: Changing default Google Maps max viewport?

texts: I want to show a portion of the world of on my Google map, where the user can zoom in, but if he zooms out, he shouldn't be able to go further out than my default viewport (i.e, if I define Africa, user should be able to zoom into Nigeria but should never see America).


texts: tutorial for deriving solar insolation estimates from high resolution Digital surface model

texts: I want to estimate solar insolation using grass gis. I have a high resolution DSM generated from Large Format Digital camera.Being a novice in this field i Would like to have a tutorial and steps to derive the solar estimates.Also I would like to estimate solar radiation on 3d building roof tops.


texts: Get attributes from many feature classes?

texts:
I have a very large .gdb (15000+ geodatabase feature classes). I know that's unusually and impractically large, but it's the result of an iterative model, which ran through every block group in NY State to calculate residential percentage based on tax parcel data. Each feature class now has one attribute -its residential percentage.

I now need a table with all of those feature classes' attributes so that I can plug the residential percentages in to an Excel document I have. Normally I would merge them and export to Excel, but of course that's not possible with such a huge .gdb.

What other way can I do this? Is there some way to derive a table from the geodatabase itself?


texts: OpenLayers: adding a loading indicator for WMS requests

texts:
I'm new to Openlayers and GIS in general. I have a WMS request followed by the addLayer method, and surrounded by the addition and deletion of a 'loading' class to an HTML element:

$(toggle).parent('div').first().addClass('loading');
Map.currentOverlay = new OpenLayers.Layer.WMS(name, Constants.webServices.wms, {
    format: 'image/png',
    tiled: false,
    transparent:true
  }, {
    singleTile:true,
    ratio:1,
    getURL: Projection.getDefaultOverlayUrl
  });
}
if (Map.currentOverlay != null) {
  Map.currentOverlay.initialZoom = map.zoom;
  Map.currentOverlay.isBaseLayer = false;
  Map.currentOverlay.visibility = true;
  Map.currentOverlay.events.on({moveend: Map.__onZoomUpdate});
  map.addLayer(Map.currentOverlay);
  $(toggle).parent('div').first().removeClass('loading');

My goal is to add a loading gif in a dialog window (with CSS) while the layer is not displayed yet. According to this post: https://stackoverflow.com/questions/9569847/when-does-openlayers-wms-physically-call-the-remote-server the request is done in the addLayer method. But it seems asynchronous: the class 'loading' is added and removed straight away, while the map is not yet displayed.
Can I add a callback or something similar to addLayer ? If not, how do I do?

NB: I didn't produce the above code, I just added the "loading" behaviour. Do not hesitate to tell me if something seems wrong with it.


texts:
I just installed QGis 2.0.1 and wanted to move the layer tools that were located on the very left hand side by default. I tried to move them to the top, next to all the other tools. Now, they just disappeared and I can't find them anywhere!

Where did they move to?

I checked all tools but still not there...

I'm very new to QGis, I already read something about Python but didn't get it and therefore didn't work. Also screen setting to lower resolution didn't work...

I would be really happy if someone could give me a very simple, step-by-step answer, where to click and what to do exactly!


texts: OpenLayers 3 XYZ layer issues with OpenRailwayMap

texts: Code with issue on Codepen - Having successfully integrated the OpenRailwayMap API into Leaflet, I wanted to use OpenLayers instead and followed the API instructions, but am getting an error: " Cannot read property 'XYZ' of undefined". Can anyone see the issue here?


texts: I'm looking for a shapefile of the Dominion Land Survey township grid for British Columbia. Where might I find this data?


texts:
I have a large dataset of multiple drill holes that have been tested. The sheet is an excel sheet and contains the holeID, depth, easting, northing and assay results.

example:

Is there a way I could get Qgis to display the highest down hole value for each drill hole?


texts: Looking for Melbourne and Canberra metropolitan suburb shapefiles / ABS coded lists

texts:
I am looking for shapefiles or a list of suburb names for Melbourne and Canberra.
I have ABS shapefiles, but I need to filter that to end up with only metropolitan suburbs.

My shapefile has codes in it for Census Geographical Area Classification Listing - I just need to know which suburbs belong to the metros.


texts:
I have a PostGIS 2.0.3 table with 113,862 records with multipolygons represented both as geometry and geography. There is a gist index on each of them. I'm running the following queries:

SELECT *, ST_AsText(geog) As wkt
FROM ds_pg_data_noise_poly WHERE
ST_CoveredBy(ST_PointFromText('POINT(12.432404 55.697712)', 4326), geog)
LIMIT 100;

The above query use Geography in a WSG84 coordinate system and returns the correct 8 rows in 33872 ms.

SELECT *, ST_AsText(geom) As wkt
FROM ds_pg_data_noise_poly WHERE
ST_CoveredBy(ST_GeomFromText('POINT(715692.377 6177776.241)'), geom)
LIMIT 100;

The above query use Geometry in a UTM32 coordinate system and returns the correct 8 rows in 10494 ms.

Is there anyway I can optimize the speed of this query? I believe I read somewhere that I could use bounding boxes to decrease the amount of records to iterate. Is this correct, and how would that look? Any suggestions are welcome.

UPDATE:

By minimizing the returned columns I managed to reduce the query time to 2168 ms.

SELECT noise_cl, noise_in
FROM ds_pg_data_noise_poly WHERE
ST_CoveredBy(ST_GeomFromText('POINT(715692.377 6177776.241)'), geom);

I still do not find this acceptable, but it is definitely better. As each record is placed within a subset, I can actually reduce the records to iterate even more.

SELECT noise_cl, noise_in
FROM ds_pg_data_noise_poly WHERE
(subset = 'subset1' OR subset = 'subset2') AND
ST_CoveredBy(ST_GeomFromText('POINT(715692.377 6177776.241)'), geom);

The above query returns in 1353 ms. I am still not quite happy with the result. Any suggestions are more than welcome.


texts: How to align points in ArcGIS Desktop 10.2?

texts:
I know this should be such a simple task and I looked around but could not find a satisfactory / dead simple answer.

I am creating a point layer and needs to draw a row of points. How would I go about making sure that these points are aligned? Something like hold shift key and it'll align it to the nearest points?

Doing this in ArcGIS Desktop 10.2.2 aka ArcMap. Layer is stored in PostgreSQL 9.2.9 / PostGIS 2.1.3.


texts: How to import geoJSON point file to ArcGIS Javascript API

texts:
I am new to ArcGIS Javascript and so far I was using ArcGIS online service to import point feature layers. But now I need to use geoJSON files to import point feature layers. How can I do that? I checked the samples in the ESRI website and could not find anything.

I am using the ArcGIS API for JavaScript v4.3.


texts: Extracting value of raster with coordinates using R?

texts:
i'm currently working on the distribution of different species on an area, for 23 years.

I'll show you a sample of my data, it will be easier to get.

My dataframe, let's call it df.

   x       y
4001758 3138416
3990685 3088576
4002641 3078682
3946723 3153793
3975356 2974350
4001284 3137528
3946723 3153793
3946723 3153793
4000195 3103181
4000168 3103446
3969985 3104761
3905824 3107504

I created a raster from .tif file.

Here is what i want to do :
I want to extract the value for each cell (with the extract() function) from my raster, where the coordinates corresponding to this value, are also in df.
I did it, so i'm still working on it, and for now, i didn't see any errors.

But the difficulty that i have is that i want to also get the coordinates corresponding for each value in my result of the extract.

Here is what i did :

test_extractcropland<-extract(cropland,df_coord)

(cropland is my raster, df_coord is a dataframe, you have a sample, above)
Do you have any idea of how i could do that ?
I looked on this forum, and stackoverflow (and others), but i didn't find something that i could adapt for my problem, or it's more likely that i didn't see how i could adapt it.

Here is some data : https://we.tl/3xe509KSZO
"coord_for_extract.csv" is what i want to filter with, the other are just raster if you want to test.

df_coord is coord_for_extract, r is my raster


texts: How to create predefined values in a field in PyQGIS?

texts:
I would like to make by programmation (python) the same thing than we can do directly in QGIS when you create a list of values for a field.
I would like to create a new field and specify a list of possible values for this field. I didn't find any function for that in the API. Is there anyone who has the solution?


texts: Georeferencing map with unknown projection and grid

texts:
It is another question about georeferencing a map with unknown projection and grid.
You can find map here.
As you can see it is already georeferenced, but poorly, I think. Projection used is ESRI:102115, Old Hawaiian UTM Zone 5N. First issue: after I gathered some info, I am not sure Zone 5N could be the right zone for Maui, rather it could be Zone 4N with ESRI:102115. So I tried to georef map with the last one, but still not a great result (error seems comparable).
Then I tried to dig more deep into it and I find some info:

Maui belongs to a series of map about Hawaii (Hawaii Division of Hydrography
Geologic and Topographic Map 1936–1960).
As you can see, in Maui there is no info about projection and grid, but I find info about projection and grid used with Hawaii. Polyconic projection, 1928 Hawaiian datum 5000 yard grid U.S. zone system, central meridian 158°.
Polyconic projection was amply used till 1950s from USGS.
Here the projection is mentioned and it is a Polyconic (this time for Maui, not for Hawaii). Same here.

I am using QGIS and I defined a new projection with proj as is:
+proj=poly +lat_0=20.54167 +lon_0=-156.3375 +x_0=0 +y_0=0 +ellps=clrk66 +towgs84=61,-285,-181,0,0,0,0 +units=m +no_defs. Then I used it to georef Maui, but still not able to achieve an acceptable result.
So I have three questions:

What is the projection used in this map?
I also need to know grid used, so I can easily georef it. What is the grid coordinates system?
Is Old_Hawaiian_UTM_Zone_5N wrong? Old_Hawaiian_UTM_Zone_4N? I think UTM is new compared to the date of the map (1942). Could be Old Hawaiian / Hawaii zone 2 (EPSG:3562) correct?


texts: ArcGIS Split Raster tool results in black rasters

texts:
I have a rather large raster that I am trying to split into small tiles in PNG format to upload to Autodesk Infrastructure.

I came across the Split Raster tool and thought it would do exactly what I wanted, so I tried using it.

The result was a series of raster in the correct place, but they were all completely black.

When I tried to add them in ArcMap they looked fine, but everywhere else they were black.

I thought maybe the raster was too large, so I tried it on a small segment and the exact same thing happened.

Anyone know what the problem is and/or how to fix this?

It would be much easier and faster than trying to figure out a way to do this programatically using clip and a series of square polygons.


texts: Data driven pages displaying same feature?

texts:
I'm creating data driven pages that display specific roads within specific counties. For whatever reason, the same roads are displayed on multiple pages even though the the correct route id is shown. I've provided pictures to help clarify. As you can see, the same route is displayed, but the route id is correct. When I zoom to the feature in data view, it shows me the correct route.

The process I used to create the data driven pages is as follows:

Add the polyline shapefile from the appropriate database connection.
Use "Select by Attributes" to identify specific roads.
Create layer from selected features.
Set up data driven pages

Add inset map, title, north arrow, scale bar, legend, and page text.

Of the 90+ counties I've done, I've only had this issue with 6 of them. Has anyone experienced this problem before?


texts: Custom coordinate 2D Helmert transformation in QGIS to Sweref 99 20 15

texts:
I have a local cartographic system in a shp file which translates to Sweref 992015, EPSG:3016 using a 2D Helmert translation. The coordinates are in the shp file x=7130.125, y=3060.129, s=0.9999555792 and theta=82284.8452 (a=0.9214369159 and b= -0.3884136588). Corresponding point in epsg3016 is x=7174622.121 and y=198905.500.
I tried using a proj string: cct -z 0 +proj=helmert +convention=coordinate_frame +x=7130.125 +y=3060.129 +s=0.9999555792 +theta=82284.8452 <Helmertin.txt >>Helmertout.txt
I am not sure if there is support still in QGIS 3.16.1 having Proj 6.3.2? But if so what have I done wrong in the proj string?


texts: How to add attributes to clustered points with Postgis?

texts:
I already have found a similar answer here:

Spatial clustering with PostGIS including attributes

but grouping after by cluster.

I cluster my data in this way:

SELECT ST_AsText(unnest(ST_ClusterWithin(geom, 2)))
from (
values
    (ST_MakePoint(10,9),'mark',20),
    (ST_MakePoint(10,9),'mark',22),
    (ST_MakePoint(10,10),'steve',21),
    (ST_MakePoint(10,13),'john',23)
) T(geom,name,age);

that correctly returns:

                        st_astext
----------------------------------------------------------
 GEOMETRYCOLLECTION(POINT(10 9),POINT(10 9),POINT(10 10))
 GEOMETRYCOLLECTION(POINT(10 13))

But what I really need is that once the data is clustered by "location" it needs to pick the first in alphabetical order desc, and if the same name, by the age desc.
A sort of:

select * from (
    select *, row_number () over (order by name desc, age desc)
    from (pre-clustered-data) T
) K where row_number = 1;

How can I do that with Postgis?


texts: Python script to read CAD data to find polygon extents

texts:
I have Civil 3D 2012 with ArcGIS for AutoCAD 300. I want to read the extents of a polygon in the drawing. Feature Class in A4A300 = Boundaries, the CAD layer is PROP-BNDY, there is only one polygon on that layer since it is the perimeter of the site.
Do I have to convert the polygon to a temporary shapefile, read it, and then delete it? I will be using the boundary extents plus a constant in all directions to clip raster images. The process should not involve user input. I have 30 directories that I will be looping through, a drawing file in each directory and a polygon in each drawing. I have the python routines for creating a list of directories that have the *map.dwg files in them and I'm familiar with clipping the rasters, I just cannot find how to read the polygon extents in the dwg file.
Your help is greatly appreciated!!


texts: How do I fix out of bounds geometric extent error?

texts:
I'm new to GIS, and I'm roughly following the instructions here:

http://mapbox.com/tilemill/docs/guides/postgis-work/

to try and import data I obtained from here:

http://factfinder2.census.gov/faces/nav/jsf/pages/searchresults.xhtml?refresh=t

into a Postgres database and then query from Tilemill. I'm getting the following error when Tilemill tries to run the query:

Detected out of bounds geographic extent (-180,4424208.8281,-13492834.8813,85.051) for layer 'acs'. Please ensure that the SRS for this layer is correct. Its native extent is '-13693889.5411,4424208.8281,-13492834.8813,4566130.4316'

Does this error have something to do with the parameter 4269 I used in the following command I used to populate Postgres?:

shp2pgsql -s 4269 140_00 census_tracts | /usr/local/bin/psql -d default


texts: Highlighting features selected by Python script in ArcMap and ArcGIS Viewer for Flex?

texts:
I have a database of orders, with geo data attached.

I want the subset of features related to a specific order to be highlighted on screen (for starters in ArcMap, but later in Flex)

in a python geoprocessing script I use SelectLayerByAttribute, and this seems to work.
(I printed the GetCount and got 11 reported).

However, I do not see the selection being hilighted in the map.
(since I have to search multiple layers later, a script seems unavoidable)

should a selection made in this way be automatically visible in arcMap (and selected in table view) ?
If not, how do I achieve this effect, i.e. that a selection made in the geoprocessing script is visible by the user firing off the script ?

code :

# -*- coding: utf-8 -*-
# Import arcpy module
import arcpy

#input parameters
Order_To_Find   = arcpy.GetParameterAsText(0)

try:
    mxd = arcpy.mapping.MapDocument(r"E:\GIS\GIS_projecten\034_Werkorders\04_Map\orders.mxd")
    whereclause = " \"ORDERNUMMER\" LIKE \'%CEM MIGRATIE_ 17%\' "
    for lyr in arcpy.mapping.ListLayers(mxd):
        print lyr.name
        if lyr.name == "MUT_E_ASSET_GEO_ACT_P":
            arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",whereclause)
except:
    print arcpy.GetMessages()


texts: UpdateCursor RuntimeError: An invalid SQL statement was used

texts:
with arcpy.da.UpdateCursor(inFeatures, "[TDS_Name]") as change_cursor:

for x in change_cursor:
        if x[0] == None:
                x[0] = x[0]
        else:
                x[0] = x[0].translate("-F", " ")

        change_cursor.updateRow(x)
    del change_cursor

RuntimeError: An invalid SQL statement was used. [SELECT [TDS_Name], OBJECTID FROM Test_Pit_2]
My purpose: Replace "-" with a " ", but only the first "-" in the ["TDS_Name"] field.
examples in TDS_Name: 1000-XX5, 4000-XX3, 3000-XX35-1
What they should be in the end: 1000 XX5, 4000 XX3, 3000 XX35-1
I used this line here which worked, but if I ran the script again would, unfortunately, remove the second "-":
arcpy.CalculateField_management(inFeatures, "TDS_Name", "!TDS_Name!.replace('-',' ', 1)", "PYTHON3", '', "TEXT")

I have imported ArcPy, and my workspace was set and the features showed up with the Exists function.
I tried using .replace in the UpdateCursor, but it didn't work and there is obviously an issue with my setup of the cursor in the first place.


texts: How do I make a single polygon in a multi-polygon feature into a new feature?

texts: I don't want to delete the polygon, but instead make it a new feature. There are over 100 polygons in the original feature, so 'split feature(s) parts' is not an option. None of the editing tools available seem to work for this particular situation.


texts: Geoexplorer - How to open query panel in a bigger window?

texts:
The query panel in Geoexplorer (displayed in lower left corner) is so narrow, that entering data in it may be cumbersome.
How can i configure Geoexplorer to open this panel  in a bigger window and with wider input (edit) control?

regards
Simon


texts: Export to CAD: How to prevent ArcGIS Pro 2.01 from adding results to Map/TOC?

texts:
How can I prevent ArcGIS Pro 2.01 from adding "Export to CAD" results to Map/TOC?
This is easily done in legacy ArcGIS Desktop in Geoprocessing options and unchecking the "Add results of geoprocessing operations to the display"

I don't see this option in Pro but perhaps there is another way this can be done?

I have 2 issues with specifically the Export to CAD tool:

Results are not grouped and when automatically added to Map/TOC this creates a mess and since I am typically not interested in adding the exported CAD data back into the map I always remove it.  Furthermore, the fact that even adding CAD data to map does not create a group layer as we are used to from legacy ArcGIS Desktop is highly inconvenient.
There appears to be a bug (have not reported this yet) - after adding CAD fields to feature class layers, then populating them values in attempts to customize the output CAD file, the Export to CAD tool succeeds but during the process of the resulting CAD layers being added to the Map/TOC the TOC crashes.


texts: Multiple Landsat 8 images shown for a single path, row, and date in GEE

texts:
I am trying to export an image above the Yenisei Delta using Google Earth Engine. I have its path, row, and date, and when I filter the Landsat 8 TOA imagecollection, I am able to find this image:

var path = 156
var row = 10;
var date = ee.Date.fromYMD(2013,7,22);

var ic = l8.filter(ee.Filter.eq('WRS_PATH', path))
  .filter(ee.Filter.eq('WRS_ROW', row))
  .filterDate(date, ee.Date(date).advance(1, 'day'));

var im = ic.first()

print(im)

var vizParams = {
  bands: ['B5', 'B4', 'B3'],
  min: 0,
  max: 0.5,
  gamma: [0.95, 1.1, 1]
};

Map.addLayer(im, vizParams)

GEE code here.

However, the map shows that there are three images (regions?) associated with this single image:

Closer inspection shows that 2. and 3. are the same image as 1, but projected incorrectly.

I am unable to export the image as it is too large (I did increase maxPixels). I exported a dozen or so others without problem.

How do I separate the true image from the false ones for export?


texts: Using Python from both OSGeo4W and ArcGIS Desktop?

texts:
This question was asked several years ago at: Can ArcGIS and OSGeo4W share the same python install?, but it seems that versions have changed since then.

I currently have ArcGIS 10.0 installed (not sure when I will upgrade to 10.1) and I would like to install OSGeo4W on my computer as well.  ArcGIS 10.0 installs Python 2.6.5 and OSGeo4W wants to install Python 2.7.2.  Can these two packages use the same Python installation despite different versions?

Otherwise, is there an easy way to install OSGeo4W with the correct Python version and keep the two Python versions completely separate on a single computer?


texts:
My code executes in three steps:

create FeatureLayer,
create QueryLayer and
join QueryLayer to FeatureLayer:

arcpy.MakeFeatureLayer_management(**{
    "in_features":features.get_dataset_fullpath(),
    "out_layer": features.get_layer_name(),
})

arcpy.MakeQueryLayer_management(**{
    "input_database":   features.get_attrib_QLD().get_workspace(),
    "out_layer_name":   features.get_attrib_QLD().get_layer_name(),
    "query"         :   features.get_attrib_QLD().get_spatial_query_definition(),
    "oid_fields"    :   features.get_attrib_QLD().get_OID_field()
})

arcpy.AddJoin_management (**{
     "in_layer_or_view" : features.get_layer_name(),
     "in_field"         : 'RDUWI',
     "join_table"       : features.get_attrib_QLD().get_layer_name(),
     "join_field"       : 'RDUWI',
     "join_type"        : 'KEEP_ALL'
})

However, while executing the step three (join), the code returns:

arcgisscripting.ExecuteError: ERROR 000229: Cannot open permits_querylayer
Failed to execute (AddJoin),

even though the layer in question has been created in step 2.

FYI: Step two is time consuming however, I thought the execution of the next statement would wait until the previous one has returned?


texts:
I'm trying to run OSGeo4W shell and the only font available is *ＭＳ ゴシック.  It's a real pain since it uses a weird "Y" instead of "\" making it hard to read commands.  Does anybody know how to activate consolas or another font?  I'm having the same problem on other command prompt programs as well.

In the registry editor I shows I have Lucida Console and Consolas as options but they don't show up in the shell.


texts: Identify button pop-up window no longer displays in ArcCatalog?

texts:
When I use the identify tool within ArcCatalog (10.2.1), I no longer get the pop-up window that provides the information about whatever feature was clicked. I'm not sure when this issue started, but at least the identify tool still does work as expected in ArcMap.

How can I get the pop-up window back in ArcCatalog?


texts: How to match a shapefile with seconds as units with other layers

texts:
I have some shapefile in coordinate seconds (longitude & latitude multiplied by 3600).  I have open Google street layer by OpenLayers plugin in QGIS.

How to match my shapefile with Google street layer in correct position?


texts:
I am trying to filter my vector layer to only show certain style layers.
This is how the first lines of my style look like:
    {
  "version": 8,
  "sprite": "../sprites/sprite",
  "glyphs": "../fonts/{fontstack}/{range}.pbf",
  "sources": {
    "esri": {
      "type": "vector",
      "url": "../../"
    }
  },
  "layers": [
    {
      "id": "someid_K/1",
      "type": "fill",
      "source": "esri",
      "source-layer": "somesourcename1",
      "filter": [
        "==",
        "_symbol",
        0
      ],
      "layout": {},
      "paint": {
        "fill-color": "#E0F2A1"
      }
    },
    {
      "id": "someid_K/0",
      "type": "line",
      "source": "esri",
      "source-layer": "somesourcename1",
      "filter": [
        "==",
        "_symbol",
        0
      ],
      "layout": {
        "line-cap": "round",
        "line-join": "round"
      },
      "paint": {
        "line-color": "#666666",
        "line-width": 1.52
      }
    }

Now I would like to add a filter to only show certain IDs. Something like this:
"filter": ["==", "id", "someid_K/1"]

But I can seem to only find filters on the layers level in the style spec.
Does such a filter exist?


texts: How to find a route given a set of noisy points

texts:
First of all, I'm a total noob to GIS problems.
Let's say I have a set of noisy GPS points (tracked while driving) along a number of streets. Is there a way to get the most likely route that generated those points (maybe using pgrouting)? In particular, sometimes, due to noise, a point lies on a crossing street and I don't want the resulting route to contain that street.


texts: How can I set the peprojection of tif file to a certain map extent and pixel size from a different file?

texts:
I have a tif file with the projection USER: 100032. And I want to reproject it to EPSG: 4326. I have previously done it in a similar way for a similar projection and used gdal.Warp(). The physical place is the same in the world (Alaska), but now it seems to change both the extent and most importantly the pixel size.
I have tried this solution, but it did not seem to work.
How can I make sure, I save the resolution, extent and the pixel size?
Here is the information:
Previously projected 100032 to 4326

CRS          EPSG:4326 - WGS 84 - Geographic
Extent      -179.7565878540000028,41.0971949020000054 : 179.9819728149999207,73.2747258849999952
Unit         degrees
Width        596
Height       130

Dimensions   X: 596 Y: 130 Bands: 1
Origin      -179.757,73.2747
Pixel Size   0.6035881890419462215,-0.2475194690999999447

New GeoTIF projected 100032 to 4326

CRS          EPSG:4326 - WGS 84 - Geographic
Extent      -179.9999232519033967,0.5107936126474186 : 179.6827057096168119,85.4634587402992878
Unit         degrees
Width        725
Height       171

Dimensions   X: 725 Y: 171 Bands: 1
Origin      -180,85.4635
Pixel Size   0.496113970981407193,-0.4967992112728179532

This is the gdal.Warp process that happens

gdalwarp -s_srs "+proj=lcc +lat_1=50.00000381469727 +lat_0=50.00000381469727 +lon_0=-107.0000076293945 +k_0=1 +x_0=0 +y_0=0 +a=6367470 +b=6367470 +units=km +no_defs" -t_srs EPSG:4326 -r near -of GTiff

I eventually need to turn these into nc files so I can calculate some things, but the latitude and longitude are for now completely off, so it will not work. I will only get 0's as the spatial resolutions are not equal.
How can I make sure the Extent and Pixel size are similar in the new projection as they were in the old projection?


texts: I have calculated mean,mode,median of all bands of an image and I want to display the bandname (B1,B2...) and their respective mean, mode,median in .csv but, I am having an issue with how to display all the information in a sequence in a table format like


texts: How to get DPI to find current scale in OpenLayers?

texts:
Is there a way to automatically retrieve the DPI value of an OpenLayers page? The goal is to automatically retrieve the map scale from this code:

function mapScale (dpi) {
    var unit = map.getView().getProjection().getUnits();
    var resolution = map.getView().getResolution();
    var inchesPerMetre = 39.37;

    return resolution * ol.proj.METERS_PER_UNIT[unit] * inchesPerMetre * dpi; }

I suppose it should be possible to consider the ScaleLine function which automatically gives the scale.


texts: Can't install Packages with Cloned ArcGIS Pro Python

texts:
I have a new installation of ArcGIS Pro 2.6 on a Windows 10 computer.  I have successfully cloned Python because I want to install third-party packages to this cloned environment and still use arcpy.
I have added this new Python installation to the Path variables in Windows and I can confirm in a Windows command prompt that the newly cloned Anaconda Python is the default Python version.

I would like to install third-party libraries.  For example, simplekml .  On a normal Python installation this would be easy with pip install simplekmlcommand. Or in this case I would think conda install simplekml.
When I try conda install simplekml I get the following error.
How do I configure my computer to allow me to install third-party packages to this cloned Python version?


texts: How to get an extent for Raster using python

texts:
I have to get an extent of a given raster so that I can save it and use the same to clip another raster with that extracted extent. As of now I am stuck in the extent part.

import arcpy
from arcpy import env

desc= arcpy.env.extent("D:\GIS @ UTD\Sem 1\GISC 6317\Lab\Lab 10\temp\Tahoe\Tahoe\Emer\erelev.grid")

print desc

The following gives an output 0 0 0 0 NaN NaN NaN NaN.

Any suggestions.


texts: Show value column in stead of key column as label in QGIS 2.14

texts:
I'm using QGIS 2.14. I have a field defined as a value relation (widget). I would
like to show the value in the linked "value column" as a label in the map in stead of the value of the "key column". Is there a way to achieve this?


texts:
We've been doing QA work on a roads dataset and I realized that many of the roads don't overlay correctly with their counterparts in Google Maps or OSM.

Currently we compare our roads with the satellite images of the area they cover and try to map them according to those images (We've been using Google Satellite View and other Satellite Views to do this). The problem with this method is that I've noticed that there is a slight shift in the position of the satellite image of the same area (I found this out by overlaying the roads on the historical images in Google Earth), I've also observed that Google roads don't always overlap with the satellite view.

Image from 2012

Image from 2015

Should I continue using the Satellite views to check the accuracy? the shift is significant enough that it will pose a problem down the years,

How else can I check the accuracy of roads and correct the issues on a large dataset?

I can't use osm or google roads to check the roads accuracy as many of the roads I have are not present in google maps or osm

I realize that the best way to do QA is to do some field testing, but that is not possible on all areas of the dataset

I am leaning towards ignoring the shift the in the satellite images altogether and just trying to maintain the relative accuracy of the roads(just trying to make sure that the roads are correctly digitized and topology is maintained)


texts:
I'm trying to combine a few Spatialite Databases into a single database, so that I can properly reference them, etc. I don't care about referencing them to one another (No foreign keys, etc), but I would like to get them all into a single database. Any suggestions on how I might accomplish this?

As an example, say I'm trying to combine the key tables from this spatialite database with the ones here. I want to ensure that the geometry information is copied.

I'm running on Windows 7, and would prefer to not use a different OS to accomplish this.


texts: Counting the number of cities within in certain polygons

texts:
I would like to count the number of smaller polygons (lakes) within three specific polygons (counties). I have been able to get the to data sets to show up in the queried table but can not get the SUM of the lakes. I believe you are suppose to use the SUM but do not know where to place it without getting an error.

SELECT c.county, a.name
FROM public.counties AS c INNER JOIN
     public.lakes AS a ON ST_Intersects(c.geom, a.geom)
WHERE c.county = 'JEFFERSON'
OR c.county = 'DENVER'
OR c.county = 'ADAMS'
GROUP BY c.county, a.name;


texts: Filling field according to value of another field with ArcPy?

texts:
I am very ill at ease with Python. I am working on ArcGIS 10.3. I am working on a table in a gdb.

I want to fill an empty field of this table with an expression (words) according to the value of a second field in the same table.

In my case, if my second field contain the word "roadway", I want to fill the empty field cells councerned by this expression "route à deux chaussées".
I know that accents are not good in programming but I need to fill my field with (french...).

My script is below:

import arcpy

fc = "D:/TESTS_TOP25/POUR_SDC/ROUTIER_TRAITEMENT/DONNEES.gdb"

#field'sname
field1 = "NATURE"
field2 = "symbology"

cursor = arcpy.UpdateCursor(fc)
for row in cursor:
    if field1 == "roadway"
        field2 == "route à deux chaussée"
    cursor.updateRow(row)


texts: QGIS 2.4.0: dissolve tool renders artifacts / errors

texts:
The dissolve tool (vector>geoprocessing tools>dissolve) renders some artifacts on shapes with high number of polygons. I’m using a shapefile with 3222 polygons and I’m trying to dissolve it using the field “municipio”, but when the dissolve tool ends calculating the new shape has some artifacts

Does anyone had a similar problem? How to work around it?

I’ve also used v.dissolve from GRASS and the dissolve polygon from SAGA, but the problem persists.
I’ve used the same shapefile on ArcGIS 9.3 with the dissolve tool. None artifact/error was produced.

Edit 1

I'm trying to use the v.clean tool as sugested by WhiteboxDev. However I'm getting this error:

  An error has occured while executing Python code:

  Traceback (most recent call last):   File
  "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\gui\ProcessingToolbox.py",
  line 181, in executeAlgorithm   File
  "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py",
  line 511, in checkBeforeOpeningParametersDialog   File
  "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py",
  line 361, in checkGrass7IsInstalled   File
  "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\ntpath.py", line 96, in join
  TypeError: object of type 'NoneType' has no len()

  Python version:
  2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

  QGIS version:
  2.4.0-Chugiak Chugiak, 8fdd08a

  Python path:
  ['C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing',
  'C:/PROGRA~1/QGISCH~1/apps/qgis/./python',
  u'C:/Users/CERCOTRICHAS/.qgis2/python',
  u'C:/Users/CERCOTRICHAS/.qgis2/python/plugins',
  'C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins',
  'C:\PROGRA~1\QGISCH~1\bin\python27.zip',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\DLLs',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\plat-win',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\lib-tk',
  'C:\PROGRA~1\QGISCH~1\bin',
  'C:\PROGRA~1\QGISCH~1\apps\Python27',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\GDAL-1.11.0-py2.7-win32.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\PIL',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\python_dateutil-2.2-py2.7.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\pytz-2014.2-py2.7.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\win32',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\win32\lib',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\Pythonwin',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win32.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\six-1.6.1-py2.7.egg',
  'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode',
  'C:\PROGRA~1\QGISCH~1\apps\qgis\python\plugins\fTools\tools',
  'C:\Program Files\SpiderOak\shell_extension_lib\shared.zip']

Not entirely sure if this error is QGIS related...
Any idea why I'm getting this error and how can I solve it?


texts: Perform Fragstats-like analyses around a defined radius of multiple points with overlapping buffers on raster file

texts:
I'm trying to do patch analyses (similar to those available on Fragstats) on NLCD landscape data (raster data). I have about 2000 points all over the state of NY.

My plan is to:
1) create a 100 meter radii around each point
2) create a new raster layer consisting of these 2000 circles.
-Everything outside the circle is irrelevant and within the circle, is the NLCD landscape data.
-Even if the circles overlap, each circle will be treated as a whole circle and a separate entity from the other circles
3) Use some patch analyses tool afterwards (ie in R, SDMTools)

I have re-projected everything on WGS 1984. Some of the data is private so I've written a code with shell data.

library(raster)
map<-raster(".tif")
data<-read.csv(".csv")
colnames(data)[1]<-"row_num"
LatLon<-data[,c("LONG","LAT")]

At this point, my relevant files are
1) map - NCLD raster
2) LatLon - coordinates .

#make LatLon into a spatial points class
llCRS<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
LL<-SpatialPoints(LatLon, proj4string = llCRS)

#perform geodesic buffer
library(geobuffer)
library(sf)
buffers_100m <- geobuffer_pts(xy = LL,
                     dist_m = c(100),
                     output = "sf")
#crop map so analyses are faster
map.sub<-crop(map2001,extent(buffers_50m))
#mask to obtain raster within these circles
rr<-mask(map.sub,buffers_100m)
plot(rr)

The result rr map that I get is a big black raster. There's no error code and I don't know what this means.

I have then tried using less points but making the radius larger so I could see what was happening.

It appears to work with one point:

buffers_test <- geobuffer_pts(xy = LL[1,],
                     dist_m = c(5000),
                     output = "sf")
map.sub<-crop(map2001,extent(buffers_test))
rr<-mask(map.sub,buffers_test)
plot(rr)

Finally, if I try using more points though,  the circles are overlapping. I can't tell if each circle is being treated independently.

buffers_test <- geobuffer_pts(xy = LL[1:4,],
                     dist_m = c(5000),
                     output = "sf")
map.sub<-crop(map2001,extent(buffers_test))
 rr<-mask(map.sub,buffers_test)
 plot(rr)

The resulting image is:

Question:
How should I address the overlapping buffer problem?

I expected an image of many circles with black around the irrelevant areas vs a giant black box. rr doesn't look outrageous :

> rr
class       : RasterLayer
dimensions  : 3226, 1202, 3877652  (nrow, ncol, ncell)
resolution  : 0.0003753802, 0.0003753802  (x, y)
extent      : -73.9758, -73.5246, 41.20613, 42.41711  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : proj2001
values      : 11, 95  (min, max)

Thoughts
-I wonder at this point whether I should just loop through each of the 2000 point (replicating the 2nd image 2000x) before performing my landscape statistics on each image individually. I don't have any great ideas and feel like I've exhausted my options.

-I've seen some similar codes online but it's not with raster data or only working with a few points.

Things that I've tried

I've posted all the steps that I've tried below because I'm not confident that this is the best method. If there is an easier method that I haven't yet tried, please let me know.

Fragstats:
Error will say it aborted but doesn't give a reason why. Aside from help manual and tutorial, most websites/youtube are in French.

ArcGIS
-I've tried zonal statistics, clip, mask, converting raster to polygon, etc. For some reason, it will show that the task is complete but the output (ie the new layer) will never appear.

Converting raster to a shapefile:
- it crashes in R (ie it keeps running forever but when I check the server usage, the CPU appears as if it's not working) and ArcGIS only produces a corner of the original raster image extent. It seems shapefiles are easier to work with and there are more tools available for doing the patch type analyses that I want.


texts: Which raster-format stores Missing, Strings, Variable Names, Multiple Bands

texts:
still looking for the ideal raster-format :-)

GeoTiff stores:
- projection/coordinate-system etc

numerical values
missing values
multiple bands

It does not allow for:

variable-names. It just does "band 1,2,3..." or "channel 1,2,3...")
string-values, only numeric

Is there a format suitable for all criteria?

THX,stn


texts: Calculating distances between consecutive points along rivers (polygon) using ArcGIS Desktop?

texts:
I have radio-tracking locations of semi-aquatic mammals. Now, I would like to calculate the distance between the consecutive locations using the river as the animals would swim (opposite to as the crow flies).

I have two different files:
a) Polypoint layer with all the tracking locations
b) Polygon layer with the rivers, streams and ponds (and some buffer around it)

These approaches seem not to work for me:
Network analyst (which seems to be most likely the easiest way): no license to create layer
The linear referencing tool to locate features on a route: My animals move once in a while a bit off the beaten water track, thus I added a buffer to the water bodies = polygon layer

I'm using ArcGIS Desktop 10.0.


texts: Problem with importing NONNA-100 bathymetry data into QGIS

texts:
Total QGIS newbie.
Tried importing Canadian Hydrographic Service Non-Navigational (NONNA-100) Bathymetric Data found here https://open.canada.ca/data/en/dataset/d3881c4c-650d-4070-bf9b-1e00aabf0a1d
Followed the Web Map Service (WMS) connection instructions found here: https://pacgis01.dfo-mpo.gc.ca/FGPPublic/NONNA_100/How_to_access_and_download_CHS_NONNA_100_on_FGP_en.pdf
Data seems to import into a QGIS but doesn't look anything like the bathymetric map.
Coordinate reference system displayed in the QGIS is: WGS 84 EPSG: 4326
QGIS version 3.14
What am I doing wrong?


texts: Is there a tool in the WebADF that will let me draw a straight line and calculate the distance between the two points.


texts:
I'm completely unfamiliar with GIS type files, I'm kind of outside the field, but I have an algorithm to test that processes pixels from such files.
In this situation, I'm looking for a way to convert a GeoTIFF file (I have several of these on disk) into one that contains no metadata, just the pixel values of all bands themselves. How can this be done?
And is there any way to extract information about how the data is stored in the file, in terms of whether it is BSQ or BI?


texts: How can I make 3D visualization of terrain data and orthophoto using GRASS NVIZ via QGIS? I am quite new with GRASS. I have only managed to get a 3D view of the terrain so far.


texts: Using erfc function in r raster package, to convert temperature data to permafrost probability?

texts:
At present, I am looking at mapping permafrost in mountains, which is modelled as a function of Mean Annaul Air Temperature. To predict statistically the probability of permafrost, you use an equation with the erfc function to predict the probablity of permafrost occuring at any given point (0-1). So far I am having no such luck with this. Converting one number through erfc is simple in Excel but I cannot seem to do it for an entire raster dataset in R, whereas this is simple to do when using simple mathematical commands (adding, multiplying etc), for the whole dataset.

Any tips?


texts: Cannot load openlayers wfs layer into Internet explorer

texts:
I am having problems loading wfs layers into IE using openlayers. My code is below. It works fine using Chrome or Firefox but fails in IE. I had a similar problem with WMS but managed to track down OL4JSFProxy which I installed and that is now working fine. However I cannot seem to get WFS working. I am guessing it is a security restriction.

vecLayerSites = new OpenLayers.Layer.Vector("Sites", {

    protocol: new OpenLayers.Protocol.WFS({
      version: "1.1.0",
      //url: "http://www.landcheck.co.uk/geoserver/wfs",
      url: "http://www.landcheck.co.uk/geoserver/OL4JSFProxy/wfs",
      //url: "http://localhost:8080/geoserver/wfs",
      // featureType: "Sites",
      featureType: "DemoSites",
      srsName: "EPSG:27700",
      featurePrefix: "stm",
      featureNS: "http://www.stmenvironmental.co.uk/stm",
    }),
    projection: 'EPSG:27700',
    strategies: [new OpenLayers.Strategy.BBOX()],
    visibility: true,
    styleMap: new OpenLayers.StyleMap({
      'default': styleContaminated
    }),
    filter: new OpenLayers.Filter.Spatial({
      type: OpenLayers.Filter.Spatial.INTERSECTS,
      value: bufferedBoundaryGeometry
    })
  });


texts: Altering aliases of ObjectID, Shape, Shape_Length, Shape_Area with ArcPy?

texts:
How can I alter the aliases of the required fields (ObjectID, Shape, Shape_Length, Shape_Area) in a GDB with Python?

For example I have this code:

rf = ur'D:\Δοκιμές\Tester'
gdb = rf+'\\Ant.gdb'
p = ['Parcels','Parts']
for a in range(2):
    RFields = OrderedDict([(arcpy.Describe(gdb+'\\'+p[a]).OIDFieldName,u'Α/Α'),(u'Shape',u'Γεωμετρία'),(arcpy.Describe(gdb+'\\'+p[a]).lengthFieldName,u'Περίμετρος'),(arcpy.Describe(gdb+'\\'+p[a]).AreaFieldName,u'Εμβαδόν')])
    for key in RFields.keys():
        arcpy.AlterField_management(gdb+'\\'+p[a],key,key,RFields[key])

But it produces an error:

Failed to execute (AlterField).

Is it possible to do what I want?

The full error:

Runtime error  Traceback (most recent call last):
  File "<string>", line 8, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 3239, in AlterField
raise e ExecuteError:
  ERROR 001600: Failed to alter field name.
  ERROR 001660: Cannot alter field name on required fields.
Failed to execute (AlterField).


texts:
My issue is to count text values from vector layer. I got a layer with streets for one of administrative region and it's names in column. It looks like this

Id   Street_name

XYZ
XYZ
ABC
XYZ . . .
ABC

And so on.

My goal is to count the frequency of each street name in this layer. So the result can look like :

Street_name  Count
ABC            3
XYZ            2
ZZY            1

I know that in ArcMap 'summarize' will be helpful. But how to do it in QgiS? Statist plugin doesn't offer such possibility.


texts: Problem with schema synchronization using Import Replica Schema (Data Management)

texts:
I want to synchronize the schema of replica. The steps I am following is documented in ArcGIS Resources:

Export the geodatabase_1 replica schema using Export Replica Schema (Data Management) tool.
Compare the schema of geodatabase_1 and geodatabase_2 replica using Compare Replica Schema (Data Management)
Import the compared schema (.xml) to geodatabase_2 replica using Import Replica Schema (Data Management)

Code sample to import is:

def import_schema(admin_sde, geodatabase_2, in_replica, schema):
    """import compared replica schema to geodatabase
    :param admin_sde: admin connection to disconnect users
    :param geodatabase_2: import schema to geodatabase
    :param in_replica: import schema for replica
    :param schema: schema xml
    """
    try:
        replica_list = arcpy.da.ListReplicas(geodatabase_2)
        replica_user = (user.name.split(".")[0] for user in replica_list
                        if user.name.split(".")[-1] == in_replica.split(".")[-1]).next()
        connected_users = [con_user for con_user in arcpy.ListUsers(admin_sde)
                           if con_user.ClientName != socket.gethostname()]
        log.info('Pausing user connection to geo-database: "{}"\n'.format(geodatabase_2))
        arcpy.AcceptConnections(admin_sde, False)
        log.info('Disconnecting users except for: "{}" and "SDE" from the geo-database: "{}"\n'.format(replica_user,
                                                                                geodatabase_2))
        for user in connected_users:
            arcpy.DisconnectUser(admin_sde, user.ID)
        log.info('Importing schema file to geo-database: "{0}" on replica: "{1}"\n'.format(geodatabase_2, in_replica))
        result = arcpy.ImportReplicaSchema_management(geodatabase_2, schema)
        return result
    except Exception as error:
        log.error(error)
    finally:
        log.info('Resuming geo-database: "{}"\n'.format(geodatabase_2))
        arcpy.AcceptConnections(admin_sde, True)

When I import the compared schema through python script, it adds the new fields to geodatabase_2 but do not delete the fields which are already been deleted in geodatabase_1. This leads to too many unwanted fields in the geodatabase_2 replica. While working on ArcMap (Distributed Geodatabase -> Import Schema Changes), it opens a Import Schema Changes Wizard where we can select the dropped fields and apply the schema changes to geodatabase_2. How to achieve this via python script? I want the same schema should reflect at geodatabase_2 replica after importing compared schema.


texts: Get a JSON return from a Overpass API call

texts:
From this call:

https://www.overpass-api.de/api/xapi?node[highway=speed_camera][bbox=-5.708215989569187,43.46669501043081,-5.605835010430813,43.588927989569186]

I'm trying to get a JSON return. Reading the API, I don't know where to put the [out:json]:

https://www.overpass-api.de/api/xapi?node[out:json][highway=speed_camera][bbox=-5.708215989569187,43.46669501043081,-5.605835010430813,43.588927989569186]


texts: How many query points fall within distance X to a set of LineString

texts:
I have a set of 10 query points and I want to know how many of them are with distance within at most 4 meters of any LineString in a given table.

Note that if a point is within 4 meter distance to multiple LineString, this should only be counted once.

I understand how I can use ST_DWithin to do the proximity checking of a query point to the LineStirng, but I am not sure as how to write an efficient SQL query to get the desired answer.

For now I can only think of 10 separate SQL statements, each either returning 0 or 1 but that is highly inefficient.

This is a clearer/cleaner version of my previous question: Preprocess OSM streets to see if query points fall "on" them


texts: Connect SQL Server to a Knack database with CARTO?

texts:
I want to have a whole database from Knack in SQL Server and for it to have scheduled updates.
Is that possible with CARTO?


texts: How to export ESRI Geodatabase format .gdb from PostGIS

texts:
I found several posts explaining how to import a .gdb file into PostGIS.
Is exporting also possible?
For example using ogr2ogr to generate a .gdb file from query result.


texts: Creating Network Flow Diagram using QGIS

texts:
I have a Telecom cable network in Geopackage format.
How do I create a Network flow diagram from which to represent the relationship between Fiber distribution hub to Network access point?


texts: ArcMap 10.3 interoperability; can't get WFS to display the right set of polygons based on my search envelope

texts:
I have a public WFS, from which I want to make a selection to save in a custom feature class.

I have made an interoperability connection, chosen the right feature classes to display, set the maximum number of features to 60 000 and set a search envelope. However, the connection doesn't display the polygons from the right area, but only the area around it (seems like). Also, despite setting a feature limit of 60 000, only 15 000 features are loaded. screenshots

As expected, Esri's online help is not a help at all.

I don't have the interoperability extension, just using what's there in stock ArcMap.


texts: ArcMap field calculator error 000989 indention error?

texts:
I'm trying to make a four character output based on concatenating a single character field (in1) with either 1 or 2 zeros, depending on the length of in2, and then concatenating that . I'm getting an indention error on this but I can't see where. I looked at the other posts with this issue and none of their solutions are working for me.

def fixRouteName( in1, in2, out ):
  if len(in1) == 2:
    out = in1 + "0" + in2
  elif len(in2)  < 2:
    out =  in1 + "00" + in2
  else:
    out = in1 + in2
  return out

__esri_field_calculator_splitter__

fixRouteName( !Exist_Partition! , !Exist_Route! , !Route! )

Here is the error message:

  Description The code block used by the Calculate Field or Calculate
  Value tool has a syntax error. This error message provided will list
  the specific Python syntax error.

  Solution This error code covers a number of Python syntax errors, for
  example, Parsing error: unindent does not match any outer indentation
  level (line 3).

  The above is a Python syntax error. Line 3 has invalid indentation.
  For specific Python issues, consult the Python documentation for more
  information, or consult the Calculate Field, Calculate Value, or
  Analyze Tools For Pro help for more information on these tools.

  When returned by a script tool or when evaluating a script tool using
  Analyze Tools For Pro, this error can also indicate that the tool's
  source script cannot be found.


texts: Where (if anywhere) can we buy a license to use Google satellite basemap as input to machine learning algorithms?


texts:
Where can I find the splash image which is shown while starting QGIS, and how to customise it?

I am developing a GIS software for myself based on QGIS, but I don't want to just start anew with QGIS's include and lib as it cost much time and energy. Instead, with the QGIS interface, I can pay more attention to create my own algorithms and add them to QGIS as plugins like sextante. Yet,I got my favoured bgimage before and still want to use it in this software I am developing.


texts: Get extents of vector shape features in PyQGIS?

texts:
i am just wondering if it is possible to iterate through vector features and return the extent per feature without clipping the vector shape?

The code looks like this and comes from the PyQgis-Documentation. How can i return the extent of every vector feature in a loop.

provider = vector.dataProvider()
            feat = QgsFeature()
            allAttrs = provider.attributeIndexes()
            provider.select(allAttrs)
            while provider.nextFeature(feat):
                geom = feat.geometry()
                # --> return extent or (xmin,ymin,xmax,ymax) ?


texts: How to copy File Geodatabase, in its entirety, while it is in use?

texts:
I'm wondering if there is any way to duplicate an in use File Geodatabase?

The Geodatabase I have is being used by ArcGIS Server for a few of our web mapping applications. I understand I could stop the services and copy the Geodatabase after, and that is an option, but I'd have to do it at 2 AM when our users aren't online.

So being lazy and wanting to enjoy my sleep, I want to find a way to copy this GeoDatabase in its entirety while it is in use.

I am a pretty good programmer, so a programming or ArcGIS solution would be perfect.


texts:
I would like to make a graph from a shapefile, where edges would express spatial relations between polygons according to DE-9IM

For that, I use fiona + shapely

My code is working but it is extremely slow because I have to check all relation for each couples

import fiona
from shapely.geometry import shape, mapping

polys = fiona.open("....shp")

def TPP(A, B):
    rel = A.relate(B)
    return rel[0] == '2' and rel[5] in '01' and rel[6] == 'F'

import itertools

# that makes 124750 couples
for pol1,pol2 in  itertools.combinations(itertools.islice(polys, 500), 2):

    # equals
    if shape(pol1['geometry']).intersects(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " equals ",\
        pol2["properties"]["ID_ILOT"]

    # touches
    if shape(pol1['geometry']).touches(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " touches ",\
        pol2["properties"]["ID_ILOT"]

    # overlaps
    if shape(pol1['geometry']).overlaps(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " overlaps ",\
        pol2["properties"]["ID_ILOT"]

    # within  a.contains(b) == b.within(a)
    if shape(pol1['geometry']).within(shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " within ",\
        pol2["properties"]["ID_ILOT"]

    # TPP
    if TPP(shape(pol1['geometry']), shape(pol2['geometry'])):
        print "poly ", pol1["properties"]["ID_ILOT"], " TPP ",\
        pol2["properties"]["ID_ILOT"]

    # disjoint
    #if shape(pol1['geometry']).disjoint(shape(pol2['geometry'])):
    #    print "poly ", pol1["properties"]["ID_ILOT"], " disjoint ",\
    #    pol2["properties"]["ID_ILOT"]

Maybe an other approach would be to use 3×3 intersection matrix with relate, which should avoid unecessary calculation.

But even with only one operation (eg intersect) it takes plenty of time for 500 polygons, and my goal is to use this method on 10000+ polygons.

Does anyone have a better approach ?

EDIT:

By using tqdm, it tooks me 10 mins to perform the calculation over 124750 couples.

100%|██████████████████████████████████| 124750/124750 [09:23<00:00, 221.26it/s]

by processing the matrix only (more or less the same duration)

100%|██████████████████████████████████| 124750/124750 [09:27<00:00, 219.93it/s]

(without any other task alongside)

If I use 10000 polygons instead, I will have 49995000 couples.
I dont think that my method is viable for such number

EDIT 2:

I figured a way to do it, I use a R-tree as spatial index to compute spatial relation on possible couple
(geopanda have a nice way to do it)
I used this tutoriel


texts:
In QGIS 3.14, I have a polygon layer that I use to show piechart - see screenshots. Everything is perfect except the drawing order of the different values (colors). I would like each pie (circle) to start at the top of the circle (noon) clockwise with the biggest value. Right now, it always has the same order: starting with orange (where applicable), than blue, red etc., as you can see e.g. in no. 14. In this case, it should start with blue at "noon", because the blue value covers almost 75% of the total amount. No 13 however is shown as I'd like to have it: orange representing the biggest class, it is also the first one.
I draw the piecharts from layer properties - charts tab, select piechart and than select the attributes I want to represent. Thus, the order is fixed. I suppose the order can be dynamically changed with the expression editor found in the tab where the attributes can be selected. However, I don't know how to write an expression that achieves what I want to do.
Thus, let's say I have the following attributes: orange, blue, red etc.
And the following features (as in the screenshot): 13, 14, 15....
The documentation is quite summary concerning diagrams:

Generated attributes with Expressions can also be used.

Unfortunately, there is no more information. As a start, I tried if ("orange" >  "blue" , "orange", "blue") - but that does not work and even if it would, defining all cases would be too cumbersome as I have 11 different attributes (colors) to be represented (however, only about 3 to 5 per pie). How can the the colors (classes) be ordered, starting with the biggest and going to the smallest value?
Edit
I tried this expression as well:
array_get (
    array_sort (
        array(
            "orange" ,
            "blue" ,
            "red" ,
            [... and so on ...] )
    ,false),
0)

This works fine to order my values acording to size, but assigns each time the same color for the first pie - bit the biggest category is not always of the same attribute (color).
Screenshot: no. 13 is correct (biggest valua at the beginning), no. 14 is not - here, blue is the biggest pie and should be the first one  - starting at noon, going clockwise:


texts: ArcGIS Mobile that uses GPS location to route to a destination?

texts: Does anyone know of an example of routing within ArcGIS Mobile 10 that uses GPS location to route to a destination?


texts: Operational layers loaded from ArcGIS Online Feature Services randomly do not display on a map

texts:
I'm developing an Android app using ArcGIS Android SDK 100.4.

I'm using a simple basemap and add feature services from ArcGIS Online to it. This has been working perfectly for a few months. Recently I had to add the same feature service urls 2 times with different symbols and this is when I started experiencing a strange behavior:
After I load the map, all layers show a LoadStatus of LOADED and isVisible() returns true. However, some of the layers are randomly not displaying on the map.
The randomness makes me think this is some race condition related to the loading order of the operational layers (= feature services), but I wasn't able to detect anything special about the layers that were not displaying.

My code for adding the layers is pretty simple:

final ArcGISMap map = new ArcGISMap(baseMapUrl);
for (String featureServiceUrl : featureServiceUrls) {
    ServiceFeatureTable serviceFeatureTable = new ServiceFeatureTable(featureServiceUrl);
    FeatureLayer featureLayer = new FeatureLayer(serviceFeatureTable);
    map.getOperationalLayers().add(featureLayer);
}
mMapView.setMap(map);

* I have omitted from my code above a part which adds a SimpleRenderer for rendering the feature services.
I'm sure I'm doing something wrong, but I'm having trouble debugging and figuring out what causes some of the layers to randomly disappear, while their status is LOADED and visible. The behavior is also not consistent (different layers are not appearing on the map every time and sometimes they all appear).


texts: Filter features by attributes from another layer

texts: I already asked this question here, but realized that my description was too meager for people to help me. Therefore, I will try to describe the problem here in more detail. Data: layer 1 (tmc_segments), layer 2(regional_roads(earth)). Layer one contains simplified geometry of the road network. Layer two contains some roads from Layer one, but these roads are more detailed. Objective: to get from layer 1 only those roads that are on layer 2. Layer one contains two distinctive attributes:road_numbe and road name. Layer two contains one attribute Name, wich contains both number and name.


texts: Clean way to install gdaltools plugin for QGIS on OSX

texts:
I have qgis 1.8.0 on osx mountain lion and am installing the gdaltools 1.2.29. I have run into a problem that a number of people are seeing.  The plugin cannot find my gdal python (see below).

Anyway, I can fix the problem hackily (if that is a word) by symlinking

  cd /Library/Python/2.7
  sudo mv site-packages/ site-packages-old
  sudo ln -s /Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site-packages /Library/Python/2.7

So that fixes the prob because the python path list that the plugin knows about includes /Library/Python/2.7/site-packages  But the real solution is to add a path to the list of paths that the python plugin[s] know about.  How do you do that?

----start of error message----

Couldn't load plugin GdalTools due an error when calling its classFactory() method

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 164, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/Users/aberezin/.qgis//python/plugins/GdalTools/init.py", line 32, in classFactory
    from GdalTools import GdalTools
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/Users/aberezin/.qgis//python/plugins/GdalTools/GdalTools.py", line 39, in
    raise ImportError( error_str )
ImportError: No module named osgeo [python-gdal]

Python version:
2.7.2 (default, Jun 20 2012, 16:23:33)
[GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)]

QGIS version:
1.8.0-Lisboa Lisboa,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/aberezin/.qgis//python', '/Users/aberezin/.qgis//python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC']


texts: Cleaning polygon after using dissolve to repair geometry in QGIS?

texts:
I'm trying to obtain the boundaries of shapefile, thus to clean the inside as a matter of speaking. While I had expected that dissolve while get me a proper polygon, unfortunately it still has some remainders. How could I get rid of them?

Does QGIS have a tool similar to Repair Geometry?

Here's the source

Result below:


texts: Recode 1 raster conditionally, based on many others

texts:
the basic premise is that I have 1 chief raster with values 55 to 75 and 6 rasters (A to F) which have all sorts of values between 1 and 15.

If the chief raster has values (e.g.) 50 to 55, those values are replaced by values from raster A in the correct and corresponding position. If the chief raster has values 56 to 58 those values are replaced by values from raster B in the correct position, etc etc.

I've gone down the path of creating a matrix for every raster then creating an index for every value range in the chief raster then subsetting 6 times.

ext <- extent(0,10,0,10)

t1 <- raster(nrows=10, ncols=10,ext)
t1[] <- sample(seq(from = 55, to = 75, by = 1), size = 100, replace = TRUE)

r1 <- raster(nrows=10, ncols=10,ext)
r1[] <- sample(seq(from = 1, to = 3, by = 1), size = 100, replace = TRUE)
m1 <- as.matrix(r1)
r2 <- raster(nrows=10, ncols=10,ext)
r2[] <- sample(seq(from = 5, to = 6, by = 1), size = 100, replace = TRUE)
m2 <- as.matrix(r2)
r3 <- raster(nrows=10, ncols=10,ext)
r3[] <- sample(seq(from = 1, to = 6, by = 1), size = 100, replace = TRUE)
m3 <- as.matrix(r3)
r4 <- raster(nrows=10, ncols=10,ext)
r4[] <- sample(seq(from = 12, to = 15, by = 1), size = 100, replace = TRUE)
m4 <- as.matrix(r4)
r5 <- raster(nrows=10, ncols=10,ext)
r5[] <- sample(seq(from = 8, to = 11, by = 1), size = 100, replace = TRUE)
m5 <- as.matrix(r5)
r6 <- raster(nrows=10, ncols=10,ext)
r6[] <- sample(seq(from = 2, to = 14, by = 1), size = 100, replace = TRUE)
m6 <- as.matrix(r6)

m <- as.matrix(t1)

idx1 <- m >=55 & m <=59
idx2 <- m >=60 & m <=62
idx3 <- m >=63 & m <=66
idx4 <- m >=67 & m <=70
idx5 <- m >=71 & m <=74
idx6 <- m ==75

m[idx1] <- m1[idx1]
m[idx2] <- m2[idx2]
m[idx3] <- m3[idx3]
m[idx4] <- m4[idx4]
m[idx5] <- m5[idx5]
m[idx6] <- m6[idx6]

t1[] <- m

Is there anything more elegant and quicker?


texts: Are there any public sources for LIDAR data for the country of Israel? Searching online has not lead me to find any sources out there.


texts: QGIS Union Error: multiple overlapping polygons created where only one should exist?

texts:
The Union function in QGIS appears to behave erroneously.

I have two shapefiles that overlap with each other, but not within themselves. They are both valid (Vector > Geometry Tools > Check Validity), and ArcGIS Union works flawlessly on both.

However, the Union function in QGIS (Vector > Geoprocessing > Union) creates a large number of erroneous duplicates. Several of them appear to have exactly the same shape. Some don't even exclude the areas where the two shapefiles overlap.

Why is this error occurring? I have not been able to find posts that describe similar problems.

Additional Oddity

When I run the Validity Check tool (Vector > Geometry Tools > Check Validity) before conducting the Union, the polygon highlighted in the example is unionized correctly. However, the problem is not resolved, as there are now other polygons that exhibit the same error.

Resources

Shapefiles are available here

Error produced with QGIS 2.18.2 on OSX 10.11.6.

Error also occurs using QGIS 2.18.2 Python bindings on Linux servers.

Update (Sept 2018): Error does not occur in QGIS 3.2

Legend for first panel: green is first shapefile (single feature highlighted in red), light blue is second shapefile, dark blue are overlaps.


texts: Spatial join if percentage of area is within

texts:
I have been performing spatial joins but there is a problem.  The geometries of both layers are slightly misaligned.  I would like to join by location where parcels are within historic_subdistricts and to keep both matching and non matching records.  The query below works but the problem is that ST_Within only works if geometry of one layer is completely contained by the other.

CREATE TABLE test_join AS
SELECT t.*, m.*
FROM parcels AS t
LEFT JOIN historic_subdistricts= AS m
ON ST_Within(t.geom, m.geom)

How would one write a query to perform a spatial join if 90% of the area of a parcel in parcels is within a subdistrict in historic_subdistricts ?

Here is a visualization of the problem:


texts:
I have a .csv-file containing x and y-coordinates as well as different variables with attributes (var1 to var3). The csv looks like that:

Loading the data into QGIS with Add Layer → Add Delimited Text Layer then results in a point layer covering the whole world with the points lying about 1/4 degree apart from each other. Filtering by variable attributes (show only points where var1 is not NULL) gives me a point layer where each point hits some parcel of land. This looks like this (point layer over some polygon layer):

What I want is to rasterize the point layer into equally sized rasters in order to be able to calculate zonal statistics and do other operations.


texts:
I have a DEM of terrain that I have successfully added a raster of woodland to:

I am trying to add to that a raster of a built up area, but the extent of the settlement raster is smaller than my terrain/woodland raster:

However, when I use raster calculator and add "terrain_woodland@1" + "settlement_raster@1", I get the two added together but my output is cropped to the extent of the settlement raster.

Can anyone tell me what I am doing wrong?


texts: Landscape software to analyze amount of dirt moved

texts:
I'm about to embark upon a design project for re-contouring landscape around my home. Assuming I get an accurate measurement of the landscape into the computer, and after I modify the landscape model to represent the proposed result, I want to produce the 'diff':

A visual representation showing the volumes where dirt needs to be removed.
A visual representation showing the volumes where dirt needs to be added.
A report of the actual volume of dirt being moved (the sums of those volumes).
A report of the delta (+/- of dirt that will be needed, or will need to be removed).

What software exists that can produce this report for me? I will prefer inexpensive software and software that is easy to learn quickly, but will accept any that meets the above criteria.


texts: Saving/exporting selected features in QGIS

texts:
I would like to save my attributes selection as a layer. I have read on other posts to do it by selecting the layer with the attribute selection, then right clicking on it and choosing "save selection as...". However, by right clicking on the layer with selected attributes, there is no "save selection as..." in my case.

How can I save/export the attribute selection as a new file?


texts: Where is conf.py file of QGIS Training Manual?

texts:
The QGIS Training Manual, at section 1.2.3, says:

"Having created your localised dataset, the final step is to replace the tokens in the conf.py file so that the appropriate names will appear in your localised version of the Training Manual."

Where and what is the conf.py file? I used "Finder" but it does not appear to be on this computer.

I am using QGIS 2.12 on Macbook Pro which runs on Mac OS X Yosemite 10.10.5.


texts: How do I insert a geometry as WKT to SQL Server Spatial with PHP?

texts:
I have a web mapping application built on Bootleaf (Leaflet + Bootstrap).  Using the leaflet Draw plugin the user can create a polygon on the map.  I'd like to insert that polygon into a geometry field to an SQL Server 2012 using PHP.  I have a field called Boundary with spatial type Geography in a table called tblProjects.

There are hard coded values here just to get things working.  They'll come from a jQuery post once I have things on the server side ironed out.

<?php
    include_once("connection.php");
    $connectionInfo = array("Database"=>"$dbname", "UID"=>$user, "PWD"=>$pass);
    $DBH = sqlsrv_connect($server, $connectionInfo);

    if( !$DBH ) {
        echo "Connection could not be established.
";
        if( ($errors = sqlsrv_errors() ) != null) {
            foreach( $errors as $error ) {
                echo "SQLSTATE: ".$error[ 'SQLSTATE']."
";
                echo "code: ".$error[ 'code']."
";
                echo "message: ".$error[ 'message']."
";
            }
        }
    }

    $sql = "INSERT INTO tblProjects (ProjectID, Boundary) VALUES (1, geography::STGeomFromText( "POLYGON (-93.30388806760311 27.994401411046173, -94.62224744260311 33.37641235124676, -79.70281384885311 31.80289258670676, -93.30388806760311 27.994401411046173)",4326));";

    $stmt = sqlsrv_query($DBH, $sql);
    if( $stmt === false ) {
        die( print_r( sqlsrv_errors(), true));
    } else {
        echo "SUCCESS!";
    }

It doesn't like my geography type, and I'm not sure how to include it; here's the PHP error:


texts: Creating geocoding/gazetteer web service?

texts:
I would like to write a geocoding web service for genealogy in Java.  Because it is for genealogy, I don't need to go down to the street level. I only need to go down to the county.

Looking at genealogy data, misspelled names are very common.  People also use a lot of non standard abbreviations.  For instance for Baltimore County might look like:

Baltimore
Baltimore County
Baltimore Co
Baltimore Cty


texts: Turfjs buffer + Leaflet return Uncaught _u {}

texts:
I'm trying to set a buffer with turfjs in a polygon. The polygon is loaded by L.geoJSON(); . When I execute bufferConcesion, this line var turfBuffer = turf.buffer(feature,value,{units:"kilometers"}); returns this error: Uncaught _u {} What might bethe problem.

This is my code:

    var featLayer = new L.GeoJSON(data, {
        style:{
            weight: 1,
            fillOpacity: 0.7,
            color: 'red'
        },
        onEachFeature: function(feature, layer) {
            layer.on('contextmenu',popUpBuffer);
        }
    });

    function popUpBuffer(e){
        var id = e.target.feature.properties.id;
        var popup = L.popup()
            .setLatLng(e.latlng)
            .setContent("

Buffer (KM): 	<input type='text' id='buffer_value'><input type='hidden' id='ca_id' value='"+id+"'/>	<button onclick='bufferConcesion()'>Go</button>

")
            .openOn(map);
    }
    function bufferConcesion(){
        var value = document.getElementById("buffer_value").value;
        if(value!==""){
            var layer = L.geoJSON(data,{
                filter:isID,
                onEachFeature: function(feature, layer) {
                    var turfBuffer = turf.buffer(feature,value,{units:"kilometers"});
                    var turfBufferLyr = L.geoJSON(turfBuffer, {style:{color:'yellow',dashArray:'5.5',fillOpacity:0}}).addTo(featLayer);
                }
            }).addTo(featLayer);
        }

    }
    function isID(feature){
        var id = document.getElementById("ca_id").value;
        if (feature.properties.id === id) return true
    }


texts: TooBig: Increasing OSRM # of coordinates

texts:
I'm trying to plan out a very specific route using OSRM.  I have 100 coordinates & now I'm running into this error.
Too many trip coordinates.  Is there a way to increase the number of coordinates I'm allowed to enter?


texts:
I'm using gdalwarp to convert a large numbers MODIS Swath (MOD11, MOD06, MOD07) products to TIF and need to mosaic the output images with same day of year.

There is a example of output data:

MOD11_L2.A2014001.1300.006.2016179212636.tif
MOD11_L2.A2014001.1305.006.2016179220648.tif
MOD11_L2.A2014002.1205.006.2016179200050.tif
MOD11_L2.A2014002.1210.006.2016179203811.tif
MOD11_L2.A2014002.1345.006.2016179203914.tif
MOD11_L2.A2014002.1350.006.2016179203951.tif
MOD11_L2.A2014003.1250.006.2016180000700.tif
MOD11_L2.A2014003.1255.006.2016180000720.tif
etc.

The images have 1, 2, 3 or 4 archives. I'm trying to use the following script without sucess:

    #!/bin/bash
    echo '==> Mosaic of MODIS Images'
    g = 001
    for f in *.tif
    do
      echo "==> Processing $f";
      gdalbuildvrt ${f/.tif/.vrt} "*A2014"+$g+"*.tif"
      gdalwarp -overwrite -of GTIFF -tps ${f/.tif/.vrt} Mosaic_${f/.tif/.vrt}.tif
      g = g+1
    done
    echo '==> Script finished with sucess'
    #done


texts:
So I just installed QGIS 1.5 using the windows installer not the osgeo4w one and now any plugins that need access to PostgreSQL databases don't seem to work, they all throw a:

"Couldn't import Python module 'psycopg2' for communication with PostgreSQL database. Without it you won't be able to run xxxxx. Please install it."

My current environment variables look like:

PATH=C:\WINDOWS\system32;
C:\WINDOWS;
C:\WINDOWS\System32\Wbem;
C:\Program Files\TortoiseSVN\bin;C:\Program Files\Mercurial\;
C:\Program Files\Console2;
C:\WINDOWS\system32\WindowsPowerShell\v1.0;
C:\Program Files\Git\bin;
C:\Program Files\Git\cmd;
C:\Program Files\Quantum GIS Tethys\apps\qgis\bin;
C:\Program Files\IronRuby 1.1\bin;
C:\Program Files\Quantum GIS Tethys\bin;

and PYTHONPATH is set to:

PYTHONPATH=C:\Program Files\Quantum GIS Tethys\apps\qgis\python

I used to have version 1.4 installed and I noticed that in the python folder for 1.4 there was a folder and file called psycopg2_psycopg.pyd, but if I take that file and put it in the python path for 1.5 it still doesn't work (didn't think it would but was worth a try)

I downloaded and tried to install psycopg2-2.2.2.win32-py2.4-pg8.4.4-release.exe which is the latest verstion of psycopg2 but because have the python that is installed with QGIS it throws a "Python 2.4 not in registry error".

Anyone know how I can get this version of QGIS running with psycopg2?


texts: Calculating distance to nearest coastline from site locations?

texts:
I want to calculate the distance to the nearest coastline from my site locations (point shapefile), and possibly the max distance as well. I am using ArcGIS 10.2 and I already have a shapefile of the coastline (lines), but I don't want to do it manually for each site using the "measuring" tool as I have hundreds of sites.

How else could I do it, perhaps using Euclidean Distance or Zonal stats?


texts: Converting UTM to Lat Long inside a Game Engine(Untiy 3d)

texts:
This is my first post on Stack Exchange.

I have exhausted my research through the game engine forums, so I am trying this forum since you guys seem to be the experts.

I have an external device that is feeding in UTM grid coordinates into the computer and using Unity 3d Game Engine: www.unity3d.com. This is a VR type project.

Im using C# and Javascript, those are the two accepted languages in Unity. I have found this code for C# online and have tried working it to my project but there seems to be some conversion into the engine errors that im struggling with. Here is the code I found, is this something I can use or will it not work in this context?

Basically I want to plug in my UTM coordinates and spit out Lat/ Long.

Hopefully it is that simple to setup and I have provided enough information about my project.

 public class NewBehaviourScript : MonoBehaviour {
     public double utmEast; //The updated UTM Easting
     public double utmNorth; //The updated UTM Northing
     public string utmZone; // The updated UTM Zone user is currently in.

     public static void ToLatLon(double utmX, double utmY, string utmZone, out double latitude, out double longitude)
     {
     bool isNorthHemisphere = utmZone.Last() >= 'N';

         var diflat = -0.00066286966871111111111111111111111111;
         var diflon = -0.0003868060578;

         var zone = int.Parse(utmZone.Remove(utmZone.Length - 1));
         var c_sa = 6378137.000000;
         var c_sb = 6356752.314245;
         var e2 = Math.Pow((Math.Pow(c_sa, 2) - Math.Pow(c_sb, 2)), 0.5) / c_sb;
         var e2cuadrada = Math.Pow(e2, 2);
         var c = Math.Pow(c_sa, 2) / c_sb;
         var x = utmX - 500000;
         var y = isNorthHemisphere ? utmY : utmY - 10000000;
         var s = ((zone * 6.0) - 183.0);
         var lat = y / (c_sa * 0.9996);
         var v = (c / Math.Pow(1 + (e2cuadrada * Math.Pow(Math.Cos(lat), 2)), 0.5)) * 0.9996;
         var a = x / v;
         var a1 = Math.Sin(2 * lat);
         var a2 = a1 * Math.Pow((Math.Cos(lat)), 2);
         var j2 = lat + (a1 / 2.0);
         var j4 = ((3 * j2) + a2) / 4.0;
         var j6 = ((5 * j4) + Math.Pow(a2 * (Math.Cos(lat)), 2)) / 3.0;
         var alfa = (3.0 / 4.0) * e2cuadrada;
         var beta = (5.0 / 3.0) * Math.Pow(alfa, 2);
         var gama = (35.0 / 27.0) * Math.Pow(alfa, 3);
         var bm = 0.9996 * c * (lat - alfa * j2 + beta * j4 - gama * j6);
         var b = (y - bm) / v;
         var epsi = ((e2cuadrada * Math.Pow(a, 2)) / 2.0) * Math.Pow((Math.Cos(lat)), 2);
         var eps = a * (1 - (epsi / 3.0));
         var nab = (b * (1 - epsi)) + lat;
         var senoheps = (Math.Exp(eps) - Math.Exp(-eps)) / 2.0;
         var delt = Math.Atan(senoheps / (Math.Cos(nab)));
         var tao = Math.Atan(Math.Cos(delt) * Math.Tan(nab));

         longitude = ((delt * (180.0 / Math.PI)) + s) + diflon;
         latitude = ((lat + (1 + e2cuadrada * Math.Pow(Math.Cos(lat), 2) - (3.0 / 2.0) * e2cuadrada * Math.Sin(lat) * Math.Cos(lat) * (tao - lat)) * (tao - lat)) * (180.0 / Math.PI)) + diflat;

         print(latitude);// Looking to display in the Console Log Lat.
         print(longitude);// Looking to display in the Console Log Long.
     }
 }


texts: My Mac always showed this notification "QGIS 3.12 (QGIS 3.10) can't be opened because Apple cannot check it for malicious software"

texts:
I am doing dialectology, and I was trying to install QGIS 3.12 (and QGIS 3.10 before) to use it in my study. Every time I opened the software my Mac always showed this notification "QGIS 3.12 ( QGIS3.10) can't be opened because Apple cannot check it for malicious software".

I also tried installing Python 3.8 hoping it would work, but it didn't and I received the same message above.

How can I install and use the software so that I can start learning the software?


texts: NULL values in QGIS attribute tables not appearing

texts:
'NULL' is not appearing in my QGIS attribute tables. Columns that used to have 'NULL' in them are now just empty.

What could be the problem?


texts: Getting values from attribute table and inserting another table using ArcPy?

texts:
I have line feature class of a river network and “FromTo” table. Lines are dissolved so they are not segmented from junctions. I need to write a Python script in ArcGIS to do following:

Begin from the first row in Hydro feature class attribute table and get its ID
Insert this ID into “From” field in “FromTo” table
Get ID(s) of line(s) which intersect(s) first line (e.g. second line
intersect first line)
Inset ID into “To” field in FromTo table

For second loop

Get ID of second row in Hydro feature class attribute table
Insert this ID into “From” field in “FromTo” table
Get IDs of lines which intersect second line (e.g. first and
third lines intersect second line)
Inset IDs into “To” field in FromTo table
So on…

Perhaps, you show a different way in order to get the same result by Python


texts: How to reduce the opacity of a Google Maps layer in CartoDB.js

texts:
I have a CartoDB.js-App where I load my vis using cartodb.createVis (I do not manually define Leaflet layers).

How can I reduce the opacity of the Google Maps base layer? I tried the following, as Google Maps has the layer index 0.

var initializeMap = function() {
    // map initialization
    cartodb.createVis('map', VIZAPI, {
            shareable: false,
            title: false,
            description: false,
            search: true,
            tiles_loader: false,
            center_lat: 47,
            center_lon: 8,
            zoom: 9,
            cartodb_logo: false
        })
        .done(function(vis, layers) {
            // layer 0 is the base layer, layer 1 is cartodb layer
            // setInteraction is disabled by default
            window.layers = layers;
            layers[1].setInteraction(true);
            // layers[1] is "all"
            layers[1].on('featureOver', function(e, pos, latlng, data) {
                updateClubInfo(data.cartodb_id);
            });
            layers[1].on('featureClick', function(e, pos, latlng, data) {
                updateClubInfo(data.cartodb_id);
            });
            layers[0].setOpacity(0.3);
            // layers[1].getSubLayer(0).setCartoCSS()
            window.overlays = vis.getOverlays();
        })

    .error(function(err) {
        console.log(err);
    });
};

This doesn't work as setOpacity is not defined on the Google Maps layer.

I also came up with the idea of creating a white, "empty" sublayer between my data layer and the Google Maps layer and set its opacity. But I don't know how to do that with the CartoDB API.


texts: How to fill a hole - with islands within?

texts:
In QGIS I have simplified a layer.
At the end of the simplification process (using GRASS v.generalize), my final layer "lost" a polygon: see white space in image:

How can I recreate a polygon which would fill that empty space?
I have tried using "Fill Ring" but that did not work.


texts: Opening a grib from the web with GDAL in Python using /vsicurl throws error on Mac but not Windows

texts:
I have a Python script that attempts to open a grib file as such:
from osgeo import gdal

url = "/vsicurl/https://ftpprd.ncep.noaa.gov/data/nccf/com/hrrr/prod/hrrr.20210505/conus/hrrr.t15z.wrfsfcf01.grib2"

gdal.Open(url)

This method has worked for me previously using the same data source. However, when I run it on my mac today I get the following error. If I run it on a windows machine, the script works.
ERROR 4: /vsicurl/https://ftpprd.ncep.noaa.gov/data/nccf/com/hrrr/prod/hrrr.20210505/conus/hrrr.t15z.wrfsfcf01.grib2
is a grib file, but no raster dataset was successfully identified.

I've tried using combinations of python 3.7, python 3.8, python 2.7, gdal 2.3.3, and gdal 3.1.4. Does anybody know how to fix this error?