texts: I did a drone flight over a forest plantation and got the DSM and DTM. Later, I was interested in calculating the height of the trees within a plot. For this, I did a subtraction between the DSM and DTM. However, when verifying the height of the trees in the CHM, these are not correct, because they vary with respect to what was reported in the field. ¿Does anyone know why this could be happening and how it can be solved?


texts:
I've a website where I'm giving the functionality of creating and editing line features. I'm using Leaflet-Draw plugin for drawing, PostGIS for storing the data and GeoJSON to display it on map.

My next task is to add the Export functionality. I've search around and found that using ogr2ogr I can convert a PostGIS table to shapefile. But this query works in FWTools, how can I run it on server. Another way I found is to use pgsql2shp. But I don't know where to run this command. I run this in Postgres console and SQL Editor but getting error.

The thing is I want to know

What is the best way in my case to export data to shapefile and how to perform it?


texts: Create an action that loads the rasters that intersect a polygon feature

texts:
Using part of the code in "Developing QGIS Action which loads a raster"

from PyQt4.QtCore import QFileInfo,QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem

s = QSettings()
oldValidation = s.value( "/Projections/defaultBehaviour" )
s.setValue( "/Projections/defaultBehaviour", "useGlobal" )
vl = QgsMapLayerRegistry.instance().mapLayersByName('catalogo_CTR')[0]
qgis.utils.iface.setActiveLayer(vl)
root = QgsProject.instance().layerTreeRoot()
group_name = "Raster layers"
group = root.findGroup(group_name)
if group == None:
    group = root.addGroup("Raster layers")
else:
    pass
fileName = '\\\\rvphnas02pw\Prodotti_Cartografici\\Ortofoto\\2012_AGEA\\TIFF_GBO\\[% A_CODICE %].tif'
fileInfo = QFileInfo(fileName)
baseName = '[% A_CODICE %]'
rlayer = QgsRasterLayer(fileName, baseName)
crs = QgsCoordinateReferenceSystem()
crs.createFromSrid(3003)
rlayer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayer(rlayer, False)
group.insertChildNode(-1, QgsLayerTreeLayer(rlayer))

s.setValue( "/Projections/defaultBehaviour", oldValidation )

How can I load the rasters that intersect a polygon feature using an action?

In my case I have a layer with the grid of the boundaries of rasters (catalogo_CTR) and a layer with administrative boundaries (Comuni).

For the "Comuni" I need an action that load the rasters corresponding to the boundaries that are intersected (i.e. 8 rasters).
The action should load raster in a group layer like the action linked at the beginning.


texts:
I am using FME's Tester to select objects that are bigger than 1 meter.

My data are in EPSG:4326.

How much is 1 meter in decimal degrees?


texts: CartoDB: Geometry code too long, Excel can't handle it

texts:
I am working on a map of regions where I want to be able to visualize some simple numbers on a periodical basis.

I have a dataset with the polygons . However those are so detailed (the geometry is incredibly long), that when I export the dataset to .csv in order to add the additional data needed in Excel (and this will be done on a regular basis), it converts to invalid code that cannot be imported.

I tried copying the geometry manually, but the limit of an Excel cell is some 35 000 characters, while one polygon in this dataset has around 150 000 characters.

What can I do to be able to store the geometry in Excel and be able to import it when I need to update the numbers I want to visualise?


texts: Time Series Extraction of Modis LAI product (every 3 days a new satellite image) in small study area

texts:
At the moment I am struggling with a time series and the extraction of MODIS LAI data. I am using the band "Lai" from MODIS/006/MCD15A3H. My study area is small, so presumably, since the pixeldata from MODIS Lai has a 500x500 spatial resolution, I only have one pixel to extract with respect to the time series of it (all 3 days one Lai value for my geometry).
I already have a similar code to extract the NDVI for my studyarea that workds, however, I don't get how to do it with the LAI data from MODIS.
Below is my code for the LAI data. It works to display it in the chart, however it seems that the scale is off, since Lai values are not higher than 8 as far as I am aware of. (NDVI would be -1 to 1) (Lai is 0 to 8ish). But in my scale it has values beyond that (up to 60), exceeding the usual range it should be displayed. I am not sure what part of the code could solve the problem.
Besides that, if I filter cloud cover (which would also smooth the line chart), there seems to be no data left. Is there a chance to filter cloud cover (for a smoother line), but still getting results?
var geometry = /* color: #d63000 */ee.Geometry.Polygon(
        [[[16.62138331919588, 48.21312568441956],
          [16.619409213360917, 48.20875013496642],
          [16.620696673688066, 48.2084927380614],
          [16.623915324505937, 48.20803514036929],
          [16.624644885357988, 48.209522317919415],
          [16.624215731915605, 48.20995130341788],
          [16.624988208111894, 48.21135263100322]]]);

//draw polygon of area interested on map//draw polygon of area interested on map
var polygon = ee.Geometry.Polygon

// which dataset + where + when
var imageCollection = ee.ImageCollection("MODIS/006/MCD15A3H")
                .filterDate('2018-01-01','2020-09-23')
                .filterBounds(geometry)
                .select(['Lai'])
                //.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));

//print collection
print(imageCollection);

//zoom to geometry
Map.centerObject(geometry, 10);

// Create palettes for display of Lai
var Lai_pal = ['e1e4b4', '999d60', '2ec409', '0a4b06'];

var Lai = imageCollection.select(['Lai']);
var Laimed = Lai.mean();

// Create a time series chart.
var plotLai = ui.Chart.image.seriesByRegion(imageCollection, geometry,ee.Reducer.mean(),
'Lai',0.1,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'Lai time series',
                hAxis: {title: 'Date'},
                vAxis: {title: 'Lai'}
});

// Display of lai values, however it seems like it doesn't do anyting afterall to display on the map
print(plotLai);

Map.addLayer(plotLai.clip(geometry), {min:0.0, max:100.0, palette: Lai_pal}, 'Lai');


texts: Storing data from WMS Server locally for ArcGIS Desktop?

texts:
I am new to GIS, and I am digitizing a map that is being provided from a WMS Server. As a result, zooming in and zooming out of the map takes loads of time due to a slow Internet connection.

I was wondering if anyone knew of any tools to download the entire file from the WMS Server to local memory, so that I do not need a fast Internet connection in order to work?

I seek links to any such tools, and steps on how to use it for this purpose.


texts:
I can't figure out why it doesn't work. At first, it worked when I'm testing using this line:

outRas= PrincipalComponents(["band2","band3","band4","band5","band6","band7"],1)

But when I tried to save the Raster layer as well, it kept returning error-like raster layer

PC= PrincipalComponents(["band2","band3","band4","band5","band6","band7"],1)
PC.save(sceneName+'_PC.tif')

I used exactly the same syntax with the test line. No error messages appeared.


texts: Choosing to use QGIS LTR or another version/release type

texts:
I am after some advice with respect to the QGIS LTR schedule.
Our IT managers have recently intervened in our use of QGIS (a State Government Authority) and are wanting to put in place a policy to only install and use the LTR release version, stating that the other releases are effectively 'Beta' versions.
Their view is from a stability and security standpoint, not from training and process mapping (as that has nothing to do with them).
I have seen and read this article - https://github.com/qgis/QGIS-Enhancement-Proposals/blob/master/QEP-4-QGIS_Long_Term_Releases.rst
Obviously the LTR has a slightly longer testing window so is less likely to have any major issues, but is their view that the other releases are 'Beta' version accurate or misled?


texts: Testing if point inside multipolygon using PostGIS?

texts:
How do I know if a point (or a coordinate lon-lat) is inside of a MULTIPOLYGON using PostGIS?

I tried with ST_Contains(), but being a MULTIPOLYGON always return False.

Also I tried with ST_Dump() for transform to a simple POLYGON but I didn't get that to work.

I used:

SELECT ST_GeomFromText(the_geom,900913)

and then:

SELECT ST_Contains("previous select result", ST_GeomFromText('POINT(lon lat)',4326))


texts: Raster Calculator Blurring Image Results

texts:
I am trying to run a NDVI analysis using raster calculator on Arg Map 10.2.1 and the results create blurry or bleeding image, almost as if the map has been dipped in water. I am using the default environment settings in the raster calculator. Both the bands being entered are very clear and have high resolution, I can see individual plants and rows in a farm field, but when I run the NDVI the resolution seems to get blurred. It almost seems like the raster calculator seems to automatically create a Guassian Blur for the image. Has anyone else encountered this and know of ways to overcome this blurring?

Also I have also noticed that when running a raster calculator, the blurring only seems to occur when subtracting or dividing image bands, but not adding or multiplying, is there a reason for this?

This is an example of the Bands resolution
An example of the bad raster output


texts: Is it possible to set the project CRS via a Python command?

texts:
Is it possible to set the project CRS of an project in QGIS via a Python command?

In my case the user should open QGIS, start my written plugin and this plugin will load spatial data which has the epsg code 31467 (Gauss Krueger Zone 3). Unfortunately QGIS default CRS is WGS84 an so I get into trouble with the units and the scale.

My first approach was to use this code:

    my_crs = core.QgsCoordinateReferenceSystem(31467, core.QgsCoordinateReferenceSystem.EpsgCrsId)
    self.iface.mapCanvas().mapRenderer().setDestinationCrs(my_crs)

This changes (after the data is already visualized) the CRS to  Gauss Krueger Zone 3. At least the espg code in the lower right corner has changed. The units are still in degrees and the scale is totally wrong. 1:5000 in Gauss Krüger is 1:528822376 in WGS 84 (in my case in QGIS). After setting the project CRS manually in the properties the project crs is really GK3 and the scale and the units are right.

So my question is, can I set the project CRS in another way than I have done it? I want to avoid the manual way.

I solved the problem by adding the following code:

qgis.utils.iface.mapCanvas().setMapUnits(0)
qgis.utils.iface.mapCanvas().refresh()

This sets the units to meter and the scale is updated automatically. The 0 stands for meters, 1 is for foot, 2 is for degree and 3 for unknown.


texts: Raster Image Tiles created using GDAL2Tiles (from GDAL 2.4.2) does not align properly with Open Street Map Tile layer

texts:
I am trying show a custom tile overlay (created from a raster image) on top of openstreetmap base layer using Leaflet.

I used one of the raster images from https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/
and used GDAL2Tiles (I have GDAL 2.4.2.) to create a tile layer.

gdal2tiles.py -p raster -z 0-8 NE1_HR_LC.tif Natural

The generated tiles are not aligned properly with openstreetmap as base layer. I have attached an image on how it appears when both of the layers appear.

Below is the code in the front end.

let minzooms =0;
let maxzooms =8;
let beginzooms =3;

var osm = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {attribution: '&copy; OpenStreetMap contributors', minZoom: minzooms, maxZoom: maxzooms});

var lyr = L.tileLayer('/tiles/{z}/{x}/{-y}.png', {tms: true, opacity: 0.7, attribution: "tiled map", minZoom: minzooms, maxZoom: maxzooms, noWrap: true });

var map = L.map('map', {
        center: [40.86667, 34.56667],
        zoom: beginzooms,
        minZoom: minzooms,
        maxZoom: maxzooms,
        layers: [osm]
    });

 var basemaps = {"OpenStreetMap": osm};
 var overlaymaps = {"Layer": lyr};

 L.control.layers(basemaps, overlaymaps, {collapsed: false}).addTo(map);

 map.fitBounds([[-180, -90], [180, 90]]);

Is gdal2tiles causing this issue?


texts: Inputting Longitude and Latitude in CSV from Geotagged, Photo Coordinates (WGS 84)

texts:
I am trying to take Geotagged picture coordinates and put them into a CSV file to display in QGIS as point values. Using Picasa, I can see that the picture coordinates are in the correct locations on the sample map. However, when I take the long and lat of the pictures and put it into a CSV file--formatted for QGIS--the points show up in incorrect locations. For example:

Latitude: 11 23\'49.9" N displayed in the photo's properties --> put as 11.23499 in CSV file

Long: 122 3\'35.9"E displayed in the photo's properties --> put as 122.3359 in CSV file

When I load it into QGIS as a layer, it shows up in a different, incorrect location. Also, I have checked my CRS to make sure it is correct. Is there anything I am doing wrong?


texts: Error #010067 when running spline interpolation

texts:
I am attempting to run a spline interpolation on a point dataset of 585 features in ArcGIS 9.3, using the spatial analyst toolbar.

The dataset holds only latitude, longitude and abundance attribute data (all as 'double' field type), in addition to the standard 'FID' and 'Shape' fields. I have successfully run IDW and natural neighbour interpolations on the data, but whenever I attempt to run a spline interpolation, Arc returns the error '#010067: Error in executing grid expression'.

The aim of this interpolation is to model fish abundance data across the North Sea and eastern English Channel, whilst producing as smooth an output as possible (hence dissatisfaction from my bosses with the natural neighbour and IDW interpolations).

I have queried the error and common responses on the ESRI forums suggest either a) illegal characters are used somewhere in the filename/attributes/file path; b) the output raster name is too long or c) the shapefile has coincident points.

I have tested these theories by firstly making a copy of the shapefile and saving it to 'C:\', with no illegal characters. This produced the same result. I also had a raster name of three characters - 10 under the 13 limit. The shapefile also has zero coincident points. Some are close, however, so I ran the spline on a 'thinned out' dataset, with only disparate points. Again, I received the error '#010067'.

Does anyone have any suggestions? I have run a spline interpolation on a test dataset that I created, just to be sure the tool wasn't in some way erroneous, but that ran without issue.


texts: Will Google Earth see my kml - on different versions?

texts:
Is there any difference in the kml script to allow different versions of Google Earth to see it ?

If I write

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.1">
<Document>

</Document></kml>

How would the kml file content change if using a different Google Earth version ?

(I would like to be sure that the kml files I create are visible on currently downloadable Google Earth versions, e.g. 5 to 7.03)

I plan to use simplekml (pyPI) to generate the files. They will use GroundOverlay, NetworkLink, Region... (as created by simplekml, which I trust to do a better job than my own brain)

I tried to find any information on version changes (for Google Earth) - but have been unsuccessful. Thank you for help.


texts: How to turn on/off all labels of all layers in QGIS

texts:
I would like to turn on / turn off all labels of any layer by one klick. Like you can do in ArcGIS according to that post:

Is there a way to turn off all labels in ArcMap without clicking on every layer?


texts: Choosing Data Type of signature file for parameter in Python Script Tool for Maximum Likelihood Classification?

texts:
I am working on a script tool for a maximum likelihood classification and I would like to set the signature file as parameter.

What data type do I have to choose in the "Add Script" dialog box in order to set a signature file as parameter?


texts: arcpy.AddJoin_management() results not writing to file in ArcMap?

texts:
This is probably just me struggling with some of the ArcGIS concepts, but I can't get the results of an AddJoin_management call to write to the layer where I've made the join.

mxd = arcpy.mapping.MapDocument('C:\\My Documents\\ArcGISmaps\\blank template.mxd')
df = arcpy.mapping.ListDataFrames(mxd)[0]

## Take geography from 'Resource library.mdb'
mapSource = os.path.join('C:\\My Documents\\ArcGISmaps\\Resource library.mdb', 'regions')

#Specify data frame & layer
bLyr = arcpy.mapping.Layer(mapSource)
arcpy.mapping.AddLayer(df, bLyr, 'AUTO_ARRANGE')

imdSource = 'C:\\My Documents\\ArcGISmaps\\IMD_02.dbf'

imdTable = arcpy.mapping.TableView(imdSource)
arcpy.mapping.AddTableView(df, imdTable)

try:
    arcpy.AddJoin_management(bLyr, "Reg_ID", imdTable, "ID", 'KEEP_ALL')
except:
    print arcpy.GetMessages()
mxd.save()

When I open the file in ArcMap 10.4, the tables have been successfully added to the TOC, but there's no evidence of the join being included. If I include getMessages() in the try: block, it suggests the join has been successful.

What do I need to do to write the join onto bLyr?

A combination of SaveToLayerFile_management and reimporting the joined table as a new layer has done the trick so I now have:

arcpy.AddJoin_management(bLyr, "ONSCD", imdTable, "MSOACD", 'KEEP_ALL')
arcpy.SaveToLayerFile_management(bLyr, "C:\\My Documents\\ArcGISmaps\\postJoin.lyr", "ABSOLUTE")

postJoin = arcpy.mapping.Layer("C:\\My Documents\\ArcGISmaps\\postJoin.lyr")
arcpy.mapping.AddLayer(df, postJoin, "AUTO_ARRANGE")


texts: Determine changes in land cover of isolated features over time?

texts:
I have a a study area of 2 Landsat scenes from 2013 which were classified for a specific natural community, and these areas have been isolated and extracted and digitized into vectors so that they can be integrated into a GIS. These communities exist as features all throughout the scene, kind of like islands, if you will. I am working with landsat data, so I want to see if these areas have grown/expanded through time, preferably over two or three decades.

What is the best way to accomplish this? I was thinking turning the natural communities polygons into a raster, using that as a mask for landsat data, and running a pixel change detection to see how they have changed. But would that give me results for if they had shrunk from previous times, since the mas would only represent current size? Or perhaps is fragstats a better option? I have no experience with fragstats. Thanks for your help and ideas

Edit for being too broad: the vectorization was for creating a detailed map for a  GIS using a basemap with much higher resolution than landsat, and the change detection is a separate part of the project. I ran a supervised support vector machine classification in ENVI to obtain the results. The images I classified were corrected for atmosphere, so I will do that with the following temporal years.

Sorry for the misunderstanding, I suppose my main question is: instead of reclassifying and performing a change detection on the whole scene, what would be the most efficient way to determine land change for only areas of interest? The areas in question are isolated swamp ecosystems, from about 1-10 hectares in size. I'd like to see if they have grown or shrunk over time. Thanks and I hope this isn't too broad or difficult to understand.


texts: Calculating percent of area changed between 2 time steps using ArcGIS Desktop

texts:
I need an image to show how many percent of area in class has been changed between 2 years

For example, in output image, red is the area that has been changed 50%, blue is 35-45%, green is <35% etc.

How can I extract this image in ArcGIS Desktop 9.2 step by step?


texts: Draw a centroid line through a line-like polygon feature?

texts:
I have all these line segments travelling between points:

I want all the one close to each other to be merged into a single feature. So I dissolved them like so:

Now I just want one line travelling through the centre of each of these line-like polygons. Any functions that do this?


texts: Why do new lines appear when transforming polygon into lines?

texts:
I'm quite new to QGIS, and I want to do something fairly simple, but I found something strange in my data.
I want to transform polygons into lines, so I can apply the buffer function afterwards to the lines (that part doesn't come into play here).
However, when I convert the polygon into lines (using Geometry Tools), new lines appear that weren't in the original map (see red square in the attached image).

Does anyone know why this could be?


texts:
Update

  Looks like it's not solely the labeling engine but something in the
  projection/reprojection, on-the-fly projection, GeoPackage or some combination that's the problem - see comment below.

The centroids produced by the labeling engine in QGIS 2.18.13 do not match up with centroids produced by the Vector-Geometry Tools-Polygon Centroids. It looks like it is using a very simple centroid algorithm, not looking at the shape but the bounding box? This leads to labels being way off from the center of the polygon. Any ideas? Is there some other setting I am missing? Nothing else is labeled. Yes, I can simply label the point output shapefile from the centroid tool, but this is not ideal. Asking here before I report this as a bug.


texts: Adding x,y coordinates from table in Google Earth Engine

texts:
This might be quite basic, but I'm struggling with adding points to a map from a CSV table of x,y coordinates in Google Earth Engine (I have over 2000 points).

I have managed to add the points as a shapefile made with QGIS, but it would be useful to know how to upload the points from a table directly.

Here is an excerpt from the "Species.csv" with coordinates (uploaded to my personal folder):

x,y
-179.7568,71.17604
-179.7104,70.95287
-179.6973,71.14408
-179.62078,70.93628
-179.43982,71.17056
-179,66.75
-173.16,64.41
-172.91,64.5
-172.187278,65.875167
-172.049306,65.785611
-171.979917,65.782306
-170.58391,63.65034
-166.783333,68.35
-166.36,65.26
-166.190167,64.982833
(...)

Any suggestions on how to import these points and add them as a layer to the map?


texts: Clip all layers of a GeoPackage (GPKG) in one step

texts: I would like to clip all layers of a GPKG-File to a given extend (canvas, mask layer, rectangle...) in one step. I could use clipping in batch mode but isn't there a nicer solution doing this in one step (SQL, QGIS)?


texts: QGIS (python) action to open folder and search for string

texts:
I am using QGIS layer property actions. Is it possible for python to open windows explorer at a directory and then search for a custom string?

Currently I can open windows explorer at the directory I want using:

import subprocess

subprocess.Popen(r'explorer /search, "D:\GIS\Infrastructure"')

But how would I search for a custom string ([%Hospital%] column in attribute file), i.e. open D:\GIS\Infrastructure and list all files that have the word "Hospital" in (as if I was to enter "Hospital" into the search bar of the folder)


texts: Is there any way to return raster tiles from a bounding box with Django raster?

texts: I've fully installed Django raster and have uploaded a GeoTIFF through the admin page. I'd like to allow the user to query the raster's tiles (stored in psql/PostGIS database) with bounding box values. Essentially I'd like only those tiles that intersect the bbox to be returned. So far I've only been able to figure out how to query with the tile index numbers, not coordinates.


texts: Configurable Map Viewer (CMV) - add tabs to widget

texts:
I m trying to add a tab panel with three tabs in a widget, but I get the following message in the debugger:

  nested divs because wipeIn()/wipeOut() doesn't work right on node
  w/padding etc. Put padding on inner div.

How can this be solved?


texts: Selecting raster cells by multiband values in QGIS

texts:
I have a raster of roads saved from a WMS, this a multiband raster.
Its red(main) an grey(back) roads.

When using the raster calculator I can select most of the red road by using: Band@1 > 200

But what expression is to be used to select the grey road?

At all grey colours, the 3 bands have the same value. My guess was that this would work, but it did not: Band@1 = Band@2 = Band@3.

I have tried a couple of other variants, but no success.


texts: Should ArcGIS Desktop users upgrade to Windows 10?

texts:
The question Will ArcGIS for Desktop run on Windows 10? was never answered, and Windows 10 is being released imminently.

Somewhat worryingly, the most recent ArcGIS Desktop system requirements page doesn't list Windows 10.

Will ArcGIS Desktop 10.x continue to work on Windows 10? Should Windows users upgrade their operating systems from versions 7/8 to version 10?


texts: trying to find max value in a field and define value as variable

texts:
I am new to scripting. I have been searching around on here and other sites, and can't seem to find something quite like what I need.

I am trying to find a max value within a field, and store that max value as a variable to get used for other processes later. I think I am defining the variable incorrectly, but I can't seem to find anything on how to define the max value of a field as a variable.

Here is what I have so far.

In this code, centerline is a variable that grabs a user defined feature class (this code is being set up as an Arc tool). SAMPINS is the name of the field that I am trying to grab the max value for. I realize that I can use statistics to get the max value (and maybe that is a better option here...), as well, but I still do not know how to define this as a variable.

The traceback error is telling me that addfields is not defined on line 179. Line 179 is the line that reads "addfields - 1 < [SAMPINS]".

addfields = 0
fieldtag = str(addfields)
with arcpy.da.SearchCursor(centerline, "SAMPINS") as cursor:
    for row in cursor:
        value = row[0]
        if value > addfields:
            addfields = value

while (addfields > 1):
    arcpy.AddField_management(centerline, "F_"+fieldtag, "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
    if addfields:
        addfeilds - 1 < [SAMPINS]
        arcpy.CalculateField_management(centerline, "F_"+fieldtag, "[SAMPINS] + [SAMPDIST] + [SAMPLF]", "VB", "")
    else:
        arcpy.CalculateField_management(centerline, "F_"+fieldtag, "0", "VB", "")
    addfields = addfields - 1


texts:
I'm searching for the nearest ten neighboring points from a given point.

When I use this query:

SELECT id, geom, ST_Distance(geography(geom), ST_GeographyFromText('SRID=4326;POINT(-0.22707 51.5444204)')) as distance
FROM cars
ORDER BY geom <-> 'SRID=4326;POINT(-0.22707 51.5444204)'::geometry
LIMIT 10;

I don't get the same result as when I use this query:

SELECT id, geom, ST_Distance(geography(geom), ST_GeographyFromText('POINT(-0.22707 51.5444204)')) as distance
FROM cars
ORDER BY distance
LIMIT 10;

Three of the 10 points returned are different. And when viewing it in QGIS it looks like the best result is the last query. I need to use index-based KNN because of performance issues.

Does anyone have any ideas?


texts:
Does anyone know which specification the PyCSW database schema follows?

From which specification comes the database fields create by PyCSW?


texts: Where can I find Ethiopian electoral constituency and district (woreda) boundaries data?

texts:
I am a PhD student in Development Studies in Oxford.

I'm trying to find a map of the Ethiopian electoral boundaries. These are not the same as the district (woreda) boundaries - districts can be in more than one electoral constituency (see this list of constituencies and the districts they cover for an example http://www.electionethiopia.org/en/amhara/167-north-wello.html).

I already have shapefiles of the district boundaries.

Any help would be amazing, even pointers to international sites that cover a number of countries and might cover Ethiopia.


texts: How to make a dynamic animated movement for the graphics on the Silverlight website?

texts:
I have some graphic points in one graphicLayer on the map and I want to move their positions to the new places every 5 seconds based on some criteria. And these movement can be visualized dynamically on the website. I'm using ESRI Silverlight API.

How should I do that? It seems the map won't update its UI until the graphics' last movement were finished. Should I create a new thread to fulfill this dynamic visualization or there is any other trick?

Thanks,

Wei


texts: How do I take a contour .shp file and create a separate layer for each elevation?

texts:
The end result is to have a .dxf file in which each different elevation resides in a different layer. If at all possible I would like to be able to do the organizing in tile mill, or through some automated process. It is possible to reorder contours into layers by hand, but that's way too slow.

Other programs I have at my disposal are:
-Inkscape
-Gimp
-Draftsite

If there's another free program that would help me I'm all ears.


texts: How do you open the Symbol Selector in ArcMap?

texts:
Normally, I think, you can open the Symbol selector by clicking any of the boxes that is shown as part of the legend below a file in Table of Contents. However, when I try this for a certain file (a tif raster file) it opens the Color Selector instead of the Symbol Selector. Why do you think this is and is there a way for me to access the Symbol Selector?

I want to access the Symbol Selector in order to put a frame around the whole raster.


texts:
I'm looking for a way to convert a classified raster into polygons based on spatial clusters within each class. For the clusters to be considered as valid I need them to consist of a minimum percentage of cells from one of the classes.

For example: An area made up of 70 % (or more) cells of class "1" will be considered as a cluster of class "1" even though the area is mixed up with 30 % cells beloning to other classes. The clustering analysis therefore should be based on the distance between cells of the same class.

Another option could be to base the clustering on a minimum number of cells within a certain class, along with a definition of a maximum search area.

For example: Within a specified area there should be 100 cells of "class 1" for it to be considered a cluster.

Most tools related to clustering seem to work only for vectors. I´ve looked at the SAGA-tool Cluster-analysis but it didn't really fit my purpose. Any ideas for solving this or which other tools that may be helpful?


texts: Writing to text file using geometry tokens of ArcPy

texts:
I am trying to write into a text file using geometry tokens but they're not being called upon.
 open(infile, "w")
 cursor = arcpy.CreateFeatureclass(fc, ["OID@", "SHAPE@", "NAME"])
  ID = 0
  for row in cursor:
      f = open("rivers_IWM.txt", "w")
      f.write(str("OID@"))
      f.write(str("SHAPE@X"))
      f.write(str("SHAPE@Y"))


texts: Sum of features from one layer to another

texts:
I have two different shapefiles and I work in Arcmap. This is what I have:

The areas in picture 1 are smaller areas and the areas in picture 2 are in an aggregated form. File 1 contains number of residents in the years 2000, 2005, 2010 and 2015. What I need is to sum the number of residents by year to know how many people are living in the larger areas. So, lets say 4 of the smaller areas make out for 1 of the bigger ones. I need to sum the number of residents from the 4 areas to know the number of residents in the 1 large one. The strange thing is; I know I've done this with the dissolve-tool before but I simply cannot remember what I'm doing wrong right now.


texts: How to add LiDAR Intensity attribute to a point feature class?

texts: I need to convert LASD (ESRI LAS Dataset) to ASCII text file. I've export the data to point feature class (LAS to Multipoint), converted multipart to single part and calculate X Y and Z.My first question is -  How do I add the LiDAR Intensity attribute to the point feature class ? My Second question Is There a better way to convert LASD to ASCII ? i tried the tool las2txt at LasTools but I still looking for an answer for my first question


texts: Output not feeding next step of ModelBuilder?

texts:
I have numerous input .txt files that I need to feed into an Arc tool (Samples>DataManagement>Features>CreateFeaturesFromTextFile) that creates Features from text files (the text file contains UTM coordinates of the vertices).

I can get this to work on a single file, but I have ~180 of them to process, and want to loop this.

I am trying to use the "Iterate Files" loop , but the output file from the Iterate Step only gives me the option of setting it as a "precondition" to feed the input of the CreateFeaturesFromText tool.


texts: Making multiple points on single location all visible at once

texts:
I am having an issue with visibility of multiple symbols at the same location. The scenario is: I am an archaeologist, dealing with groups of objects all buried together at the same time ("hoards"), which of course all have the same spatial location. Hoards can be composed of objects of different materials, e.g. gold, bronze, amber objects etc.
I am trying to make a map of distribution of different materials, and I have used categorized symbols to indicate different materials, but only the uppermost symbol is visible. I thought of possibly making the symbols different sizes so that a bit of each symbol would peek out, but that implies quantitative significance. I thought maybe there might be a way to slightly offset each symbol, or maybe to make a tiny pie chart? Any help?


texts: Accessing selected features in the sublayer output from tracing geometric network

texts:
I have a geometric network set-up with stream lines and water users. I am trying to use ArcPy to automate tracing downstream of each node/junction in my geometric network to get userIDs for all downstream nodes. I can't figure out how to access the selected nodes in the layer group that is output from the tracing. Arcpy.searchcursor("sublayer name") does not seem to be able to find my sublayer.  How can I access just the selected nodes from my tracing?

When I use ListLayers as follows, it doesn't include the trace output layer group with selections even if I can see it in the TOC:

    arcpy.TraceGeometricNetwork_management(myNetwork,outNet,mypts,"TRACE_UPSTREAM")
    layers = arcpy.mapping.ListLayers(mxd)
    for layer in layers:
       if layer.isGroupLayer:
          if layer.name == arcpy.mapping.Layer(outNet):
             for subLayer in layer:
                if subLayer.name == userLayer:
                #Get a list of downstream user IDs
                   downstreamUsers = []
                   rows = arcpy.SearchCursor(userLayer)
                   for row in rows:
                      downstreamUsers.append(row.HYDROID)
                   del rows, row
           print "downstream user IDs", downstreamUsers

So I also tried referring specifically to the trace output layer and then the SearchCursor says it can't find my sublayer with the user nodes:

  arcpy.TraceGeometricNetwork_management(myNetwork,outNet,mypts,"TRACE_UPSTREAM")
    layer = arcpy.mapping.Layer(outNet)
    for sublayer in layer:
       if subLayer.name == userLayer:
          #Get a list of downstream user IDs
          downstreamUsers = []
          rows = arcpy.SearchCursor(userLayer)
          for row in rows:
             downstreamUsers.append(row.HYDROID)
             del rows, row
      print "downstream user IDs", downstreamUsers


texts: Deprecation warning when importing Processing in pyQGIS 3

texts:
I use QGIS 3.4 with Python 3.7 and I have a warning message when importing Processing like this :

import processing
from processing.core.Processing import Processing

Here is the message :

Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working

What is the new way to import it without getting this DeprecationWarning ?


texts: Making CSV used as query for selecting some features?

texts:
layer = qgis.utils.iface.activeLayer()
selected_features=layer.selectedFeatures()
print layer.name()
l = []
for i in selected_features:
    l.append(i.attribute("ID"))
print(l)

When I run this in Python console I got [10,11,12,14] now I want to export this text [10,11,12,14] in one line of csv-file, and use it later for selecting features.

With a result of print(l) I want to select features in my activeLayer(). How can I achieve this?


texts: Using Python to discriminate Personal Geodatabases and MS Access Databases

texts:
I need to list a server's geodatabases, to analyze and process them. My users are however creating both personal geodatabases and MS Access geodatabases.

Listing mdb files reports both types of databases, I woudl like to narrow the list down to geodatabases only.
The boolean workspace property is_geodatabase would be exactly what I need to use, as it indicates whether a database is spatially enabled or not. Unfortunately, this property is only available for SDE databases.

I am using ArcGIS 10.0. and the script below. The choice was made to list mdb files by extension using os.walk() rather than arcpy.ListWorkspaces() for speed of process, and because ListWorkspaces retrieves MS Access databases anyway.
How can I amend my script to report only personal geodatabases, and not flat MS acess databases?

for path, subdirs, files in os.walk(analyzed_root_folder):
    for x in files:
        if x.lower().endswith(".mdb") == True:
            ext_mdb.append( os.path.join(path, x))


texts: Drawing circle of 3 km radius in ArcGIS API for JavaScript?

texts:
I am able to draw the circle but I want a circle of radius 3 kms

I have checked all the previous posts.

I want to know what is the units of radius in the below code.

var radius, circle, ring, pts, angle;
var currentPoint = new esri.geometry.Point(longitude, latitude, this.parent.esriMap.spatialReference);

currentPoint = esri.geometry.geographicToWebMercator(currentPoint);
circle = new esri.geometry.Polygon(this.parent.esriMap.spatialReference);

ring = []; // point that make up the circle
pts = 40; // number of points on the circle
angle = 360 / pts; // used to compute points on the circle

for (var i = 0; i <= 360; i = i + 10) {
    var radian = i * (Math.PI / 180.0);
    var x1= radius * Math.cos(radian);
    var y1= radius * Math.sin(radian);
    var circlePoint = new esri.geometry.Point(x1, y1, this.parent.esriMap.spatialReference);
    circlePoint = esri.geometry.geographicToWebMercator(circlePoint);

    points[i] = new esri.geometry.Point(
        parseFloat(circlePoint.x) + parseFloat(currentPoint.x), parseFloat(circlePoint.y) + parseFloat(currentPoint.y),
        this.parent.esriMap.spatialReference);
    ring.push(points[i]);
}
circle.addRing(ring);
//rest code to add the layers to map


texts: Understanding use of spatial indexes with RTree?

texts:
I'm having trouble understanding the use of spatial indexes with RTree.

Example: I have 300 buffered points, and I need to know each buffer's intersection area with a polygon shapefile. The polygon shapefile has >20,000 polygons. It was suggested I use spatial indices to speed up the process.

SO... If I create a spatial index for my polygon shapefile, will it be "attached" to the file in some way, or will the index stand alone? That is, after creating it can I just run my intersection function on the polygon file and get faster results? Will intersection "see" that there are spatial indices and know what to do? Or, do I need to run it on the index, then relate those results back to my original polygon file via FIDs or some such?

The RTree documentation is not helping me very much (probably because I'm just learning programming). They show how to create an index by reading in manually created points, and then querying it against other manually created points, which returns ids that are contained within the window. Makes sense. But, they don't explain how that would relate back to some original file that the index would have come from.

I'm thinking it must go something like this:

Pull bboxes for each polygon feature from my polygon shapefile and place them in a spatial index, giving them an id that is the same as their id in the shapefile.
Query that index to get the ids that intersect.
Then re-run my intersection on only the features in my original shapefile that were identified by querying my index (not sure how I'd do this last part).

Do I have the right idea? Am I missing anything?

Right now I'm trying to get this code to work on one point shapefile that contains only one point feature, and one polygon shapefile that contains >20,000 polygon features.

I'm importing the shapefiles using Fiona, adding the spatial index using RTree, and trying to do the intersection using Shapely.

My test code looks like this:

#point shapefile representing location of desired focal statistic
traps = fiona.open('single_pt_speed_test.shp', 'r')

#polygon shapefile representing land cover of interest
gl = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('class3_aa.shp', 'r')])

#search area
areaKM2 = 20

#create empty spatial index
idx = index.Index()

#set initial search radius for buffer
areaM2 = areaKM2 * 1000000
r = (math.sqrt(areaM2/math.pi))

#create spatial index from gl
for i, shape in enumerate(gl):
    idx.insert(i, shape.bounds)

#query index for ids that intersect with buffer (will eventually have multiple points)
for point in traps:
        pt_buffer = shape(point['geometry']).buffer(r)
        intersect_ids = pt_buffer.intersection(idx)

But I keep getting TypeError: 'Polygon' object is not callable


texts:
I am using ArcMap 10.4.1. I am trying to create a custom color ramp that only has 4 colors, as defined by the company I work for. I am trying to create the color ramp through the Style Manager but the pre-set colour ramp has 13 spaces for color however I only need 4.

Does anyone know a simple way I can create and save a custom colour ramp with only 4 colors?


texts: Calculate intersection of line from known origin with DTED

texts:
I'm working on calculating a point on the ground where a camera is pointing from a drone.
I've got DTED2, and can extract the position / altitude of the drone as well as the az/el angles of the camera boresight.
What I'd like to return (ideally using python) is a lat/Lon and elevation for the point where the camera boresight intersects the DTED. Very little experience of working with DTED so be gentle!


texts:
Consider the following scenario:

I have made polygon of mainland Africa + Madagascar.
From a different source I have got some point data about this area.
The point data include points in small islands that I want to discard. So I use over() in the R-package sp to filter out such
points.
However, hundreds of points at the coastline is slightly off my polygon(s), so they would also be discarded.

In the image below, only the points that are off are plotted.

I think that if I could only enlarge the two polygons (mainland Africa + Madagascar) slightly, then I could filter out points using over(). Perhaps there is another solution that gives the same result, and if so that is great.

EDIT: using rgeos::gBuffer() as suggested by @mdsumner, the problem is now solved, which can be seen in this graph, where the points the are truly off mainland Africa and Madagascar are colored red and the points that where "false positives" are colored black.

The code used to enlarge the polygon(s) was:

Africa.2 <- spTransform(x=Africa, CRSobj=CRS("+init=epsg:32662"))
library(rgeos)
Africa.larger <- gBuffer(Africa.2, width = 50000)


texts: Setting rendering order of line layers in QGIS

texts:
I have a dataset where different road types are in separate layers, i.e. motorways, primary roads and secondary roads. In turn, each layer is subdivided by embedded codes into different sub-types - 4 lane roads, 2 lane roads etc.

I currently have the layer style set so that each road is outlined in black. Thanks to another answer to a question here on GIS SE, I already know how to make different road segments 'blend together' so the black outline remains 'underneath' the colour using symbol levels, and that works fine.

However, because for example, the primary roads are in an entirely separate layer to the secondary roads, and are obviously above them in the rendering order, the black outline of the primary roads cuts over the colour line of the secondary roads. This has the effect of making any secondary road that meets a primary road look like a dead end.

The drawing order is essentially:

Primary road colour
Primary road outline
Secondary road colour
Secondary road outline

When it needs to be:

Primary road colour
Secondary road colour
Primary road outline
Secondary road outline

So that the colour line of the Secondary road blends seamlessly with that of the Primary road. Within a layer this is obviously handled by Symbol Levels, but this is overridden by the layer order. Drawing Layer 0 in the Primary road layer is still higher than Drawing Layer 1 in the Secondary road layer.

How can I fix this?


texts: Exporting multiple, consistently sized raster files

texts:
I'm working with a computer that uses multiple ASCII raster files and ran into an issue with the file alignments. As you can see in ArcMap, the rasters are being aligned correctly.

However, upon exporting them to ASCII and loading them into the software they are slightly askew and this can be verified when looking at the headers,

$ head -n 6 x.asc
ncols         7752
nrows         4384
xllcorner     -90.415008381648
yllcorner     45.952043272731
cellsize      0.000348592492
NODATA_value  -9999
$ head -n 6 y.asc
ncols         8581
nrows         5190
xllcorner     -90.519812091
yllcorner     45.799266992357
cellsize      0.000348592492
NODATA_value  -9999
$ head -n 6 z.asc
ncols         8581
nrows         5190
xllcorner     -90.519812091
yllcorner     45.799266992357
cellsize      0.000348592492
NODATA_value  -9999

I've checked all of the typical problems (e.g., projections, clipping, etc.) and everything is being projected to NAD 1983 for the program. The library we are using for processing is fairly basic so the alignment needs to come from ArcGIS which it should be up to. There is also data present to the boundaries for the region so the theory was that should have been enough.

Perhaps there is a command I'm forgetting about, but how can I ensure that the raster files are all aligned correctly upon export?


texts: Some labels disappear depending on zoom (scale dependent vis is turned off)

texts:
for one layer I have an issue with a very few select labels disappearing depending on zoom level. See label "3" on the right in below screenshots.
"Scale dependency" is inactive and "Show all labels" is active. This only happens with a few labels in this layer - and always the same. I have another layer that is very similar but based on a different variable and there all labels reliably get rendered.

I found that some labels keep appearing/disappearing just by panning the map around in the same zoom level. I've played around with the various placement options and found that 'placing in whole centroid' vs 'placing in visible centroid' leads to those faulty labels not appearing at all in any zoom or panning state.

Is there a fault with the polygons/centroids?

In another layer using the same map/polygons this problem does not occur.


texts:
I found something strange in how GeoServer serves WMS, specifically how it draws donut polygons.  When the whole polygon is inside WMS GetMap request BBOX, the inner hole of donut polygon is filled. When part of the polygon falls outside WMS request BBOX, the inner hole is not filled.

Is it a bug, a feature or am I doing something wrong?

Tested on GeoServer versions 2.3.3, 2.8.3, 2.9.1.


texts: How do I make a slippy map with a GeoTiff?

texts: It seems that most applications I use for serving map data, gdal2tiles or Geoserver takes the supplied GeoTiffs and chops it up into png files with no metadata.  What is the standard practice to place these tiles correctly.  Do I just lose elevation and bounding box data and that's it?  Do I just place the tiles based off of the transform of the previous zoom level?  I feel like there is something fundamental I'm missing.


texts: Missing option to apply symbology by Categories in ArcGIS

texts:
I have a layer file with several thousand records and a lot of columns. I'm trying to map this by breaking out values from one of the columns into categories. But I don't see an option for Categories in the symbology tab in "Layer Properties". I don't know why this is the case since there are definitely values in multiple columns that could be broken into categories. Here is what the Symbology tab looks like:

Is ArcGIS being weird or am I doing something wrong here?


texts: Loading Geoserver WMTS layer in ArcGIS Desktop?

texts:
I am using GeoServer 2.13. I am trying to load WMTS layer into ArcGIS Desktop, but WMTS layer is not loading into it. When I am trying to load into QGIS, its loading successfully. WMTS URL is below

http://localhost/geoserver/gwc/service/wmts?

I also checked the User access log, ArcGIS requested URL is different than QGIS requested URL.

ArcGIS request log:

  0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:18:37 +0530] "GET
  /geoserver/gwc/rest/wmts/nurc:Pk50095/raster/EPSG:4326/EPSG:4326:13/1965/8789?format=image/jpeg
  HTTP/1.1" 401 1253 0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:18:37 +0530]
  "GET
  /geoserver/gwc/rest/wmts/nurc:Pk50095/raster/EPSG:4326/EPSG:4326:13/1964/8788?format=image/jpeg
  HTTP/1.1" 401 1254 0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:18:37 +0530]
  "GET
  /geoserver/gwc/rest/wmts/nurc:Pk50095/raster/EPSG:4326/EPSG:4326:13/1964/8789?format=image/jpeg
  HTTP/1.1" 401 1254 0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:18:41 +0530]
  "GET
  /geoserver/gwc/rest/wmts/nurc:Pk50095/raster/EPSG:4326/EPSG:4326:13/1967/8789?format=image/jpeg
  HTTP/1.1" 401 1171

Qgis request log:

  0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:30:11 +0530] "GET
  /geoserver/gwc/service/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=nurc:Pk50095&STYLE=raster&FORMAT=image/jpeg&TILEMATRIXSET=EPSG:4326&TILEMATRIX=EPSG:4326:12&TILEROW=983&TILECOL=4391
  HTTP/1.1" 200 926 0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:30:11 +0530]
  "GET
  /geoserver/gwc/service/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=nurc:Pk50095&STYLE=raster&FORMAT=image/jpeg&TILEMATRIXSET=EPSG:4326&TILEMATRIX=EPSG:4326:12&TILEROW=984&TILECOL=4397
  HTTP/1.1" 200 7429 0:0:0:0:0:0:0:1 - - [21/Jul/2018:13:30:11 +0530]
  "GET
  /geoserver/gwc/service/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=nurc:Pk50095&STYLE=raster&FORMAT=image/jpeg&TILEMATRIXSET=EPSG:4326&TILEMATRIX=EPSG:4326:12&TILEROW=982&TILECOL=4398
  HTTP/1.1" 200 7281

Can anyone suggest me how to load WMTS layer into ArcGIS Desktop?


texts: Converting text field to date field in ArcGIS Pro

texts:
How to convert date, which is in text formated field to date formated field?
Also, there is dates which have only year, without days and months and empty fields.
I've tried Convert Time Field. But it gave me null values


texts: Relationship between sp package in R; MapWindow; ArcGIS/QGIS?

texts:
I am using ArcGIS at office for some visualization work, but I feel it kinda dump to performing the same tasks several times for visualizing some data on different maps. I am think of automating it using R.

I have come across several things, sp package for R, ArcGIS for visualization, MapWindow for drawing maps, and another opensource Quantum GIS. But I puzzled on the relationship between. For example, if I use MapWindow to draw the maps, can I use only R for visualization work? And what's difference between R+sp package and Quantum GIS + python + R?

I am trying to visualize the location of beds in an elderly home which the elderly is infected with some bacteria, the bed will be shown red if the patient gets the bug, and shown green if the patient does not.

It is perfectly OK for me to use only ArcGIS if I only need to visualize the clustering of bacterial infection in one elderly home. But the problem is: I have 33 of them. And I don't like repetitive tasks to be done manually.

Can I use R plus some kind of GIS system to do the repetitive task, and give me the 33 charts?


texts: Intersection of multi-polygons in PostGIS


texts: FME convert Esri Line Shape to Esri Point Shape

texts:
How to convert an Esri Line-Shape to an Esri Point-Shape using FME Desktop?
The distance between the resulting points should be 1 meter.
Which transformer(s) would be best?
Edit

The distance between the points is not everywhere exactly 1 meter.


texts: Link to an individual info window in CartoDB?

texts:
I just published a very large parcel dataset on CartoDB.  Another website would like to link their parcel dataset to ours.  Meaning, if you click on a parcel on their site and want to know more information about that parcel that's published on our site, they can link directly to the individual parcel.  That info is shown in a CartoDB info window on our site.

Is there a way to link directly to that individual pop-up -- perhaps in the URL?


texts:
I have a Model in ModelBuilder, and I did not do any special manipulation, but when I am running Model I am getting error:

  All processes have already been run

I tried to remove intermediate data and re-run ArcMap, but it did not help.


texts: GRASS r.slope.aspect via QGIS 2.18 vs QGIS 3.4 - not working

texts:
I seem to be unable to replicate actions with GRASS tools when I access them through QGIS 3.4. In particular the r.slope tool. I receive error messages and it fails in QGIS 3.4.
They seem to work fine when accessed via QGIS 2.18.
Is QGIS 3.4 known for such issues?
I run Windows 10.

***2018-12-18T00:08:30     CRITICAL    Traceback (most recent call last):
          File "C:/PROGRA~1/QGIS 3.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 413, in processAlgorithm
          Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
          File "C:/PROGRA~1/QGIS 3.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 411, in executeGrass
          startupinfo=si if isWindows() else None
          File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\subprocess.py", line 756, in __init__
          restore_signals, start_new_session)
          File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\subprocess.py", line 1155, in _execute_child
          startupinfo)
          OSError: [WinError 87] The parameter is incorrect***


texts: Accidentally deleted URL to Story Map in ArcGIS Online - How to recover?

texts: While editing the story map's properties in ArcGIS Online, I mistakenly modified the story map's URL and now the app cannot open. Is there a way to retrieve its correct URL? For reference, this was modified under My Content -> selected story map -> Edit -> Properties section


texts:
I'd like to do a classification using  a multispectral raster (1 Blue 475, 2 Green 560, 3 Red 668, 4 Red Edge 717, 5 Near IR 840). I'd like to classify the crops (kiwi plants).

I tried a few unsupervised classification with ERDAS like NDVI, DVI, SARVI, TBDVI RDVI. I used RDVI to to do a 90 class classification, then blackened the classes that didn't show the crops.

Is there any better way to get the crops classified?


texts: Getting zone from UTM in order to convert to latlong

texts: I am trying to convert UTM to a latlong point. To get that done I need to know the zone where the UTM coordinates are. how can I find the zone from UTM?. the place coordinates that I work on lay on zone 39 and 40.


texts: How to query and display layer using ArcGIS Python script tool?

texts:
I have database which contains point shapefiles and raster files for 3 species.
My aim is to develop a script tool that should have a search tab to enter the species name.
also, a check box to select what layers to display.
how to develop this.

give me hint to develop python script for querying and display a layer.
also, i should add the script tool as a button in the toolbar.


texts:
I couldn't find any proper example of how to import a .csv with WKT coordinates (Polygon) in PyQGIS.

Here are more details on my .csv. I am working with polygons. The delimiter is a semi colum, and the name of the variable that contains the relevant coordinates is WKT . Thus I wrote :

1) uri = "data/idf.csv?delimiter=%s&crs=ESPG:2154&wktField=%s" % (";", "WKT")

2) idf_stats_csv = QgsVectorLayer(uri,"idf","delimitedtext")

But this does not work. Do you have any more advice?


texts: How to randomly replace a given line vector multiple times within a polygon?

texts:
I have this line vector that represents the movement trajectory of a given animal. I would like to randomly replace this vector line within a polygon, multiple times, changing its position and orientation, but without changing the trajectory itself. Is this possible? It is OK if random trajectories intersect each other, as long as they don't fall outside the polygon. Here is an example:

The black square is the polygon, the continuous blue line is the real trajectory, and the three dotted blue lines are examples of three random placements of the trajectory within the polygon.

How can I do this on QGIS or R?

I don't have any experience with Python.


texts: Comparing yearly changes in TaxLot land use in NYC (mapPLUTO) - tax lots change through time

texts:
I am trying to build a database for NYC in which, for specific boroughs of the city, I track changes per tax lot unit in terms of land use and other attributes, both qualitative and quantitative. For example, I am trying to track the changes in land use or Floor Area Ratio per tax lot through time.

For this, NYC offers extensive data in their database mapPLUTO, which can be found here, both for 2018 and historic data that traces back to 2002.

So my initial approach is to arrange the attribute tables of each year and compare the fields I am interested in, from year to year, creating something like a series of columns that express whether there has been any changes between year a and a-1.

This would be an easy and straightforward process if tax lots where always the same, but between 2002 and 2018 construction and development has taken place, increasing the number of lots by subdividing them. Below, two images compare the tax lots of 2002 and 2018 over an area that was developed:

This situation gets even worse regarding the coding of the data: in 2002, the lot was identifiable via 2 attributes, being block (169 in this case) and lot (the lot number within the block, being = 1). However, in 2018, all the new tax lots get a new lot number, and one of the new entities has the value 1, meaning that the datasets from 2002 and 2018 have both an entity identified as 169-1, but this entity is not the same.

My question is, then, whether there might be a methodological solution for how to be able to track changes between datasets that have different entities in place, which makes it hard to compare entity-to-entity.

I have a possible path to overcome this problem, which you may comment on, or add an alternative. My first option is to perform the analysis per block. That way, the first step of my analysis would be to ask whether the block has the same amount of tax lots or not. Then, for each block, data is retrieved from its tax-lots and aggregated at the block level, which becomes the research unit.

Another way of fixing this issue might be to do some kind of spatial join in which the latest version of the database "inherits" the info of previous years. however, the two datasets do not fully overlap even though they have the same projection, leading to some lots not to overlap with their past selves (see image below) representing 2002 and 2018, with the latter in yellow lines).

This study is being carried out in R.


texts: Changing opacity of map outside feature in OpenLayers 3?

texts:
I want to change opacity of map outside feature (which is closed polygon) in OpenLayers 3.

Is it possible?


texts: Find out whether geography point lies within geography bounds

texts:
I have three points:

Point A
Point SW
Point NE

I want to know if Point A resides within the bounds determined by points SW (for south-west coordinates) and NE (for north-east coordinates).

I've so far tried both ST_Intersects and && to determine this, which works fine for small regions, but fails whenever the bounding area determined by SW and NE points covers more than half the globe. In this case it seems to instead use a different, smaller, bounding area using the same four points.

As an example, if I'd want to cover everything from the SW point of South Africa to the NE point of Ireland it would instead catch points where the the NW point is in Ireland and the the SE point is in Madagascar. That is, China would not considered to be within the bounds.

Here are some SQL queries to illustrate this:

select
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))') as expectTrue_1,
ST_Intersects(ST_GeographyFromText('SRID=4326;POINT(0 0)'), ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))')) as expectTrue_2,
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  180 80, -80 80, -80 -80, 180 -80))') as expectTrue_3,
ST_Intersects(ST_GeographyFromText('SRID=4326;POINT(0 0)'), ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  180 80, -80 80, -80 -80, 180 -80))')) as expectTrue_5,
ST_GeographyFromText('SRID=4326;POINT(-90 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))') as expectFalse_1,
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((1 -1,  -1 -1, -1 1, 1 1, 1 -1))') as expectTrue_5;

I'd want for this to return t t t t f t but it returns f f f f t t.

I've also tried ST_Intersects(ST_MakeBox2D()) as per the following SO post:
https://stackoverflow.com/questions/24212355/query-by-coordinates-takes-too-long-options-to-optimize/24218861?iemail=1&noredirect=1#24218861

I'm using PostGIS 2.1.3 with the Geography data type.

What functions or approaches should I use to have PostGIS treat my more-than-half-globe-spanning polygons as more-than-half-globe-spanning ones instead of 'optimizing' them to be less?

I understand I could just split the bounds in, say, two parts, and check for each, but there must be a cleaner solution.


texts:
I'm using Postgis 2.1.5 on Ubuntu 15.04. In my system there is a second disk mounted with the raster data as /media/name1/data/rasters/*.tif.

The *.tif are good looking when I watch the metadata with gdalinfo.

The data is imported into the raster data table with raster2pgsql without problems.

Variables for raster support are set in the PostGIS enviroment file correctly.

When I want to clip the raster data by:

SELECT ST_Clip(a.rast,b.geom)
FROM public.raster a, public.umringe b
WHERE a.filename = b.datei_bez;

The following error is the result :

  ERROR:  rt_band_load_offline_data: Cannot open offline raster: /media/ausbild/2ndDrive/Geodaten/Komserv/Daten/LK_WIL/Work/test/raster_new/BPlan.07231001.2.0.plan.tif
  CONTEXT:  PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN
  ********** Fehler **********

  ERROR: rt_band_load_offline_data: Cannot open offline raster: /media/ausbild/2ndDrive/Geodaten/Komserv/Daten/LK_WIL/Work/test/raster_new/BPlan.07231001.2.0.plan.tif
  SQL Status:XX000
  Kontext:PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN

I tested it with the referenced file on the "normal" home-path, and then it works. But for the work with the whole 45GB of raster-data the volume is not big enough, so I have to use the "external" media drive.


texts: Issues in projection while defining data frame extent by coordinates

texts:
I'm trying to set the following coordinates to a dataframe:

custom_extent = {1:-87.678956, 2: -87.611652, 3:41.867129, 4:41.954932}

where 1 is key for XMin, 2 is key for XMax, 3 is key for YMin, and 4 is key for YMax

However, when I try to set the new extent up, the coordinates are not the results I expect. Here's a code snippet:

        print data_frame.extent.XMin                                        # prints 480718.855003
        original_spatial_reference = data_frame.spatialReference
        newExtent = data_frame.extent
        data_frame.spatialReference = arcpy.SpatialReference("WGS 1984")
        print data_frame.extent.XMin                                        # prints -87.6986700746
        newExtent.XMin, newExtent.XMax = custom_extent[1], custom_extent[2]
        newExtent.YMin, newExtent.YMax = custom_extent[3], custom_extent[4]
        data_frame.extent = newExtent
        print data_frame.extent.XMin                                        # prints -93.500018925
        data_frame.spatialReference = original_spatial_reference
        print data_frame.extent.XMin                                        # prints -2.38749449402

Every printout does not match the desired XMin: -87.678956. Any thoughts?


texts: Unscramble RGB raster - 4D interpolation/lookup

texts:
I have a 4-band raster image (RGB), downloaded from a WMS, representing clay content of topsoil ranging from brown (high clay content) to yellow (low clay content). The 4th band is alpha-band, with constant 255 values.

Instead of the red, green, blue values (all of them ranging from 0 to 255) the concentration is needed (ranging from 0 to 100). As a help, at certain points, I have the corresponding percentage (see the image above). That means, I have a lookup table like this:

R,    G,    B ,   percentage

244,  229,  115,  9

242,  187,   62,  26

220,  148,   37,  33

183,   97,    6,  41

177,   87,    1,  43

...

The question is, how can I get the value of clay content in percentage for every raster pixel.

For me this looks much like an interpolation task; however, resample, interp, approx work only in fewer dimensions.

Here is a similar problem however with a different color ramp.


texts: Clipping raster by feature using GDAL in QGIS?

texts:
I have successfully selected the feature from a shapefile (boundary) and intend to use this as a mask to clip the raster map underneath. However, having followed instructions when I open the GDAL > Raster extraction > Clip raster by mask layer I don't have an output layer and only an input. I have tried several options but without success.

Error code is:

QGIS version: 3.8.0-Zanzibar
Qt version: 5.11.2
GDAL version: 2.4.1
GEOS version: 3.7.1-CAPI-1.11.1 27a5e771
PROJ version: Rel. 5.2.0, September 15th, 2018
Processing algorithm…
Algorithm 'Clip raster by mask layer' starting…
Input parameters:
{ 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'INPUT' : '/Users/MartynMac/My Documents/GIS Projects/Tithe & Enclosure Maps/OS Georef/Ludham Catfield Crop_modified.tif', 'KEEP_RESOLUTION' : False, 'MASK' : QgsProcessingFeatureSourceDefinition('parish_region_6a39b3b8_cbb3_4593_b1de_4f956d31c796', True), 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : None, 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None }

GDAL command:
gdalwarp -of GTiff -cutline /var/folders/cf/_4glv9gn5171yhgxywymh9ww0000gn/T/processing_4c8576f1d012446297b0c24783abf44b/8e3fc98821ff46718579079472e8e09c/MASK.shp -crop_to_cutline "/Users/MartynMac/My Documents/GIS Projects/Tithe & Enclosure Maps/OS Georef/Ludham Catfield Crop_modified.tif" /var/folders/cf/_4glv9gn5171yhgxywymh9ww0000gn/T/processing_4c8576f1d012446297b0c24783abf44b/42b858586ae74639b5c42371a2a19e25/OUTPUT.tif
GDAL command output:
/bin/sh: gdalwarp: command not found

Execution completed in 0.19 seconds
Results:
{'OUTPUT': '/var/folders/cf/_4glv9gn5171yhgxywymh9ww0000gn/T/processing_4c8576f1d012446297b0c24783abf44b/42b858586ae74639b5c42371a2a19e25/OUTPUT.tif'}

Loading resulting layers
The following layers were not correctly generated./var/folders/cf/_4glv9gn5171yhgxywymh9ww0000gn/T/processing_4c8576f1d012446297b0c24783abf44b/42b858586ae74639b5c42371a2a19e25/OUTPUT.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts:
Currently, I have a strong knowledge ArcGIS/Python/JavaScript and basic RDBMS. I'm currently taking some database theory/sql classes and trying to figure out what course I should take next, I'm considering Java.

I've been told ArcObjects and .NET is a thing of the past and to focus on Java, is this true? I see most GIS/Developer job postings requesting proficiency in ArcObjects .NET C#.


texts: Programmatically "locking" map layout in composer of QGIS 3?

texts:
I am programming a plugin for Qgis 3 from which I want to export a composer layout containing different maps that are supposed to have different extends and show different layers.

I am having trouble "fixating" a map before I make the changes required for the next one. When I use the composer interface, there is the option to use "Lock layers", which seems to do what I want. By looking into template files saved from the composer and the API documentation, I have figured out that QgsLayoutItemMap.setKeepLayerSet() is probably the corresponding method in python. The only problem is that it simply does not seem to work. The maps are not "fixated" and everytime I change something in the map canvas, it is still reflected in all maps contained in the layout.

Does anybody have any ideas what the problem might be?

    #Accessing the main map by recasting it as a QgsLayoutItemMap
    main_map = sip.cast(composition.itemById("Main"),QgsLayoutItemMap)
    main_map.zoomToExtent(raster_layer.extent())
    main_map.setKeepLayerSet(True)

    QgsLayoutExporter(composition).exportToImage( output_folder + "/test1.png", QgsLayoutExporter.ImageExportSettings())
    #Here the exported map does not show the background layers

    #Activating the background map
    root = QgsProject.instance().layerTreeRoot()
    for c in root.children():
        if c.name() == background_map.name():
            c.setItemVisibilityChecked(True)

    QgsLayoutExporter(composition).exportToImage( output_folder + "/test2.png", QgsLayoutExporter.ImageExportSettings())
    #And here it does...


texts: Updating Arc Toolbox paramaters with validation script

texts:
I'm simply trying to disable parameter 1 if parameter 0 if entered. I modified the following, but it does not work, do I need to return anything?

def updateParameters(self):
  """Modify the values and properties of parameters before internal
  validation is performed.  This method is called whenever a parmater
  has been changed."""

  if self.params[0].value == False:
      self.params[1].enabled = False
  else:
      self.params[1].enabled = True

  return [self.params[0],self.params[1]]


texts: Different values between original raster and cropped raster with rasterio.mask

texts:
I have a large sentinel-2 image and I want to clip it according polygons which are stored in a geojson file.

My code is:

df = gpd.read_file(os.path.join(geojson_folder,geojson_file))
df.head()

    id  name    geometry
0   1   almond  (POLYGON ((749159.7552032885 3643270.487528285...
1   2   olive   (POLYGON ((749639.4218857342 3643450.362534203...
2   3   banana  (POLYGON ((749055.617041968 3643819.579651611,...
3   4   orange  (POLYGON ((748500.2135149258 3643882.693688774...
4   5   apple   (POLYGON ((747988.9898138981 3643119.013839092...

Where gpd is geopandas, and the crs is espg:32636

To perform clip (using rasterio.mask) for each polygon by each sentinel-2 raster (format jp2) I use:

os.chdir(output)

coords = [json.loads(df.to_json())['features'][i]['geometry'] for i in range((df.shape[0]))]
field_name = [crop for crop in df.name.tolist()]
for poly,fn in zip(coords,field_name):
    os.mkdir(fn)
    for sn2_folder in os.listdir(src):
        if sn2_folder.endswith('.zip'):
            pass
        else:
            for sn2_band in os.listdir(os.path.join(src,sn2_folder)):

                data = rasterio.open(os.path.join(src,sn2_folder,sn2_band))
                epsg_code = int(data.crs.data['init'][5:])
                srs = osr.SpatialReference()
                srs.ImportFromEPSG(epsg_code)

                try:
                    out_img, out_transform = mask(dataset=data, shapes=[poly], crop=True, all_touched=False)
                    out_meta = data.meta.copy()
                    new_name = sn2_band
                    out_tif = os.path.join(fn,new_name)
                    out_meta.update({"driver": "JP2OpenJPEG",
                         "height": out_img.shape[1],
                         "width": out_img.shape[2],
                         "transform": out_transform,
                         "crs": srs.ExportToProj4()}
                                 )

                    with rasterio.open(out_tif, "w", **out_meta) as dest:
                        dest.write(out_img)
                except:
                    pass

print ("finish")

The outcome seems satisfying. I get clipped rasters, which are correctly geo-referenced, with the desired pixels size etc. However, I've noticed that the pixel values of the original raster are different from the clipped raster. Why is that happening? is there some hidden resample (I didn't see anything in the official documents)?

UPDATE:
After playing with this some more, the (weird) solution is to change the "JP2OpenJPEG" to "Gtiff". For some reason this creates jp2 files with matching values between original and clipped rasters.


texts:
Using geoserver 2.9.2. I'm generating vector-tiles (geojson output, EPSG:27700) for the full extent of a layer, for only one zoom level. Estimated number of tiles is around 17.000.000. After full generation, when I can see that each job has completed, let's say, 2125000 tiles, on disk I have only around 1.000.000 files.

the disk quota is set at 10 GiB, and the total size of tiles cache is 4.5 GiB.

I tried to use only one thread... same result: the tiles completed displayed by GWC interface is different by the number of files on disk.

In OL, when I hit the missing tile, and on the request header I can see the MISS, on disk, the missing file is generated.

The geoserver is running in a docker container, with cachedir mapped on a NAS.
(I tried also on a fresh installed geoserver (not using a docker container)

Somebody know what is with this difference between GWC completed tiles and actual files on disk?

note1:
if for a given extent is generating 16 PNG files, for the same extent is generating only 1 geojson file.

The ratio between PNG generated files and GEOJSON generated files is 16:1,
which makes sense: for ~17.000.000 estimated tiles to generate ~1.000.000 files.
The rest of vector tiles are missing, and when I request a missed one, the header is marked as MISS and is generated on disk.


texts:
I recently changed my OS to Xubuntu 16.04.
Now i try to install GDAL to use it via shell and QGIS.

I followed the steps from the official UbuntuGIS website (https://launchpad.net/~ubuntugis/+archive/ubuntu/ppa). While updating I can see the path for the UbuntuGIS ppa (I post the log ant the end).

Now my Questions is: Can I use UbuntuGIS under Xubuntu and if yes, what am I doing wrong?

Output of the shell while updating my ppa:

Hit:1 http://ppa.launchpad.net/ubuntu-desktop/ubuntu-make/ubuntu xenial InRelease
Hit:2 http://de.archive.ubuntu.com/ubuntu xenial InRelease
Get:3 http://security.ubuntu.com/ubuntu xenial-security InRelease [102 kB]
Hit:4 http://de.archive.ubuntu.com/ubuntu xenial-updates InRelease
Hit:5 http://ppa.launchpad.net/ubuntugis/ppa/ubuntu xenial InRelease
Hit:6 http://de.archive.ubuntu.com/ubuntu xenial-backports InRelease


texts:
I have a plpgsql function that has to return the result in a new table. Below is the code

CREATE OR REPLACE FUNCTION vrg.VRGTest_Ver5(_tbl1_VRG regclass,_tbl2_urban regclass,_tbl3_postcodes regclass,t_name varchar(30))
RETURNS SETOF record AS

      $BODY$
      BEGIN

      EXECUTE format(

     ' CREATE TABLE IF NOT EXISTS %I.%I
(
    PostcodeName text,
    vrgid double precision,
    urban_gid integer,
    vrg_area double precision,
    urban_area double precision,
    weights double precision
 )',
    'vrg',
    ('' || t_name)
    );

       EXECUTE format  ('SELECT
             postcodes.name as PostcodeName,
             fishnet.fid as vrgid,
             urban.gid as urban_gid,
             st_area(fishnet.geom) as vrg_area,
             ST_AREA(ST_INTERSECTION(urban.geom, fishnet.geom)) as urban_area,
             ST_AREA(ST_INTERSECTION(urban.geom, fishnet.geom))/st_area(fishnet.geom) as weights
         FROM '   || _tbl1_VRG       ||  ' AS fishnet,'
              || _tbl2_urban     || ' as urban,'
              || _tbl3_postcodes || ' as postcodes
         WHERE
             ST_INTERSECTS(urban.geom, fishnet.geom) and ST_INTERSECTS(fishnet.geom, postcodes.geom)')
         into t_name;
       END;
       $BODY$ LANGUAGE PLPGSQL;

How can I modify the code to write both create  and select statement in execute.format or is there any better way? I'm new to plpgsql.
When I run the function

  select vrg.VRGTest_Ver5( 'vrg.bb_sumatrapart_grid_1000m' ,'vrg.ghs_pop_indonesia_v2_reclass_poly_union_v2_projected','vrg.gfk_2016_id_5_digit_pcd_areas2013_projected',
                          'vrg.ResultTest')

it goes to infinite loop.


texts: Adding new feature class to already classified image in ArcGIS Desktop?

texts:
I have at hand this already classified image with 8 vegetation classes.  Now I would want to change a part of existing vegetation class to a NEW vegetation class, making it 9 classes.

I do not want to reclassify the image as it will take some time.  Also, I just needed to change few polygons for the 9th class.

How can I do this?

Is there a way to change the class manually?

I have ENVI 5.2 and ArcMap 10.2.

This is a portion of the original image with 8 classes:

This is a portion of the new 9th class:

This is how I want it to be (In here, as just did an overlay for visualization purposes. Can the CON tool add the new class 9 to the original image?)

When I applied CON, the output was this, which was just the same as the 9th class.

Finally, with help from @MichaelMiles-Stimson, the merging works.  A new class 9 was added to an already classified image with 8 classes by doing this in the Raster Calculator: Con(IsNull("rip_raster"),"classif9.dat",9)


texts:
I have an atlas in QGIS 2.14 with a filtered legend whose position is anchored at the bottom right corner. When viewing atlas pages individually within the print composer, the legend size increases with number of items on each map, but does not decrease when the number of legend items is reduced on subsequent pages. On export, the legend box is at the same maximum size on all pages, leaving a gap under the item list where there are less visible layers. I don't remember having this problem in previous versions - is this a bug in 2.14?

Edit: I have logged this as a bug - issue 14707.


texts: Service Area - Not taking into account of all network

texts:
I'm trying to work out an 800m service area from a station.
However, QGIS is not taking into account some of the network lines.
The yellow colour highlights the network that should be included in the service area (service area highlighted in orange). Anyone got any idea about what could be wrong?  Please note network lines not included in the services area is the proposed road network.

Network Data - Modified Open Street Map Data shp from Geofabrik. Modified to include the proposed road using Civils 3D and converting it back into shp. I have tried to draw directly into QGIS the missing road network lines, however, it still comes out with the same result.

I have redrawn the lines and crossing them into existing road so it is clearly intersecting but the results are the same.


texts:
I am working on a river watershed on QGIS and I need a land cover map.

I have a good raster of land use data with 15 m resolution for year 2015, but it covers around the 90 % of the basin.

I have also global land cover raster maps with a resolution of around 300 meters for year 2012 and 2013.

Is there a way to generate a raster data covering the missing 10% of the 15 m raster by using some geostatistics on the 15 m available data and the 300 m raster?

In the case I would assume the land use of 2015 equal to the 2012 and 2013.


texts:
After upgrading to 7.0.4, I could not execute r.mapcalc perfectly, there is a problem :

syntax error, unexpected INTEGER, expecting VARNAME or NAME
Parse error
ERROR: parse error
ERROR: An error occurred while running r.mapcalc

It happened when i try to mosaic some images, all the required parameters were already set up correctly.

Anyone have a clue ?


texts: Runtime error arcpyy using fields[f].append(row.getValue(f))

texts:
I am trying to write a arcpy script that will turn off all NULL and empty fields.  The script runs for a while and then gives a run time error at fields[f].append(row.getValue(f)).  Below is the script.

import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True
arcpy.env.addOutputsToMap = False
arcpy.env.workspace = "C:/temp"
out_layer = "temp.lyr"

mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, '*')[0] # need to find way to get list of layers there can be more than one data frame
listfields = []
thebadlist = []
for LayerNeedsFieldsTurnedOff in arcpy.mapping.ListLayers(mxd):
    if LayerNeedsFieldsTurnedOff.isGroupLayer:
        print "group print"
    elif LayerNeedsFieldsTurnedOff.isFeatureLayer:
        fields = dict((f.name, []) for f in arcpy.ListFields(LayerNeedsFieldsTurnedOff) if not f.required)

        rows = arcpy.SearchCursor(LayerNeedsFieldsTurnedOff,"","","","")
        for row in rows:
            for f in fields.keys():
                fields[f].append(row.getValue(f))  #**the error happens here**
                #print row.getValue(f)

        for field, values in fields.iteritems():
            #print field

            if  field == "Shape":
                listfields.append(field)
            elif all(map(lambda s: s is None or not str(s).strip(), values)):
                thebadlist.append(field)
            else:
                listfields.append(field)

        # LayerNeedsFieldsTurnedOff = arcpy.mapping.ListLayers(mxd, 'exampledata', df)[0]
        # fill in your desired fields to remain visible

        desiredFields = listfields

        field_info = arcpy.Describe(LayerNeedsFieldsTurnedOff).fieldInfo # I need this for field count

        for i in range(field_info.count):
            if field_info.getfieldname(i) not in desiredFields:
                if field_info.getfieldname(i) == "SHAPE":
                    field_info.setvisible(i, 'VISIBLE')
                elif field_info.getfieldname(i) == "Shape":
                    field_info.setvisible(i, 'VISIBLE')
                else:
                    field_info.setvisible(i, 'HIDDEN')

        arcpy.SaveToLayerFile_management(LayerNeedsFieldsTurnedOff, out_layer)
        arcpy.MakeFeatureLayer_management(LayerNeedsFieldsTurnedOff, 'temp_layer', '', '', field_info)
        refLyr = arcpy.mapping.Layer('temp_layer')
        # rename the ref layer the same as your target layer
        refLyr.name = LayerNeedsFieldsTurnedOff.name
        lyrfile = "C:/temp/temp.lyr"
        print LayerNeedsFieldsTurnedOff.name
        arcpy.ApplySymbologyFromLayer_management(refLyr, lyrfile)
        arcpy.mapping.UpdateLayer(df, LayerNeedsFieldsTurnedOff, refLyr, False)
        # mxd.save()
        print 'cleaning up-'
        if arcpy.Exists('temp_layer'):
            print '\'temp_layer\' still in memory...deleting now...'
        arcpy.Delete_management('temp_layer')
        print 'deleting obj refs...'
        del LayerNeedsFieldsTurnedOff, refLyr,  fields
        print 'done.'


texts:
I have read similar threads on here regarding the technicality and am in way over my head. I am not a professional cartographer - I am a professional coach driver that uses google maps to plan routes.
What I would like to know is why, when you scale out, countries seem way out of proportion. E.g. The UK and Ireland is the size of half of Europe.


texts: Smoothing borders of polygons taking a layer of points into account

texts:
I'm trying to smooth borders that separate clusters of points.
GRASS instrument v.generalize.smooth works pretty well, except that in some cases the points end up on the wrong side of the smoothed border. On the picture: green line is the correct "unsmoothed" border, purple line -- is the "smoothed" border; the bottom point should belong to a different cluster than the upper three points, yet, due to smoothing it is now in the same cluster.

My question: is there an instrument in QGIS (or anywhere else) that allows to smooth borders of polygons taking a layer of points into account?


texts: Change CRS identification in file in raster without reprojecting QGIS

texts:
I have a series of rasters saved as geotiff's that are assigned the incorrect CRS.  Is it possible to change the CRS of the raster without reprojecting?
When I open the files, QGIS reads the CRS information from the file as EPSG 2285 (Washington State Pane North, US Feet), and puts them in the wrong place.  The rasters are actually in EPSG 2286 (Washington State Plane South US Feet). Changing the assigned CRS to EPSG 2286 moves the rasters to the correct location.
I would like to save the correct CRS information with the file.  Is it possible to do this without reprojecting?  If it were a shapefile, I could just change the .prj file so it matches how the file is actually saved, is something similar with a .tiff raster?


texts: Accounting for spatial autocorrelation in rasters

texts:
Background:

I have a set of fine resolution rasters for the Florida peninsula. These rasters are going to be used in species distribution modeling. Many of these rasters represent variables that are clearly going to be related, so prior to using them in the model I'm performing a global PCA to generate new rasters. All of my rasters are of continuous data (no binary data or categorical variables).

However, I think I can make the case that it doesn't make sense to perform a global PCA for such a large area where the values of the variables are going to vary dramatically from one part of the state to another. So, I was considering going with a geographically weighted PCA (gwpca) instead, and was thinking that it might make sense to base the radius of the gwpca on the distance of spatial autocorrelation. Ideally, I would determine the distance individually for each pixel, but I don't know of any software that can do this, any paper or website that outlines the process, or even if it's computationally viable if I were to try and code my own solution. Alternatively, the backup plan is to come up with a single estimate and use it across the raster. This spatial autocorrelation distance would also be used to thin out high concentrations of points in my point data that is going to be used for generating the model.

Question:

How do I determine the distance at which spatial autocorrelation occurs? Should I do it on a per pixel, per site, and/or per raster basis?

What I've done so far:

Most of my work is done in R (I also use tools like GDAL and GRASS. I'm on Linux, so cannot use ESRI gis products).

I can take a raster and find the global Moran's I. I can also subset it using any given point and a radius, and use this to generate semivariograms and correlograms. I have a basic understanding of sill, nugget, and range, but am not aware of any how I could automate the process of determining the distance in which spatial autocorrelation is a factor.


texts: Using a color override for representations

texts: I'm looking for a way to automatically assign a color to a feature, mainly polylines.  I want to do this by populating a field with a value so I can link it to a color.  I believe this is possible with representation symbology, but I haven't been able to find out how to use the override field for color.  Is there a specific type of field it has to be?  Is there some format the attributes have to be, ie RGB values?  Perhaps there are other ways to do this as well?  I am open to ideas, but representation symbology is what I need the end result to be.


texts: How to classify entire pixels within loaded polygons as a feature collection, in Google Earth Engine

texts:
I need to determine the Modis pixels that enter whole within the polygons contained within a feature collection.

I also need to determine those who don't fully enter.

Apologies for English

This is the script

https://code.earthengine.google.com/2c58d1f0e127f0d002ed20cf6872fe05


texts: I have offset some data points in QGIS but I want my points to stay aligned in the legend? It looks funny when points aren't lined up. Anyway I can do this?


texts: Determining distance from Source in ArcMap?

texts:
I am trying to determine the distance in kilometers of a large number of points from a polygon which is their source. The polygon indicates a geological quarry and the points (which number more than 900) contain stone material from this quarry. I am aware that I can manually determine the distance with the ruler tool, but that would be time consuming considering there are over 900 points.There must be a tool which can determine the distance from the source for each point at the same time. Ideally what I am looking for is an attribute table that displays the kilometer distance for each point.I am using ArcMap 10.3.1. I

am fairly new to this software.


texts: arcpy add rasters using raster calculator

texts:
what I want to do is add multiple raster layers stored in FGDB together using the raster calculator in arcpy

To start off I have seen this question Using Loop with Raster Calculator in ArcPy?

this is what it would look like in arcgis raster calculator

'remap1'+'remap2'+'remap3'+'remap4'+'remap5'+'remap6'+'remap7'+'remap8'+'remap9'+'remap10'+'remap11'+'remap12'+'remap13'+'remap14'+'remap15'+'remap16'+'remap17'+'remap18'+'remap19'+'remap20'+'remap21'+'remap22'+'remap23'+'remap24'+'remap25'+'remap26'+'remap27'+'remap28'+'remap29'+'remap30'+'remap31'

now in Python

import arcpy
from arcpy.sa import *
arcpy.env.workspace=env
reclass1 = '+'.join(["Raster("+"r'{}'".format(x)+")" for x in arcpy.ListRasters() if x[:5]=='remap' if int(x[5:]) < 14 ])
print reclass1

which prints out to

 Raster(r'remap1')+Raster(r'remap2')+Raster(r'remap3')+Raster(r'remap4')+Raster(r'remap5')+Raster(r'remap6')+Raster(r'remap7')+Raster(r'remap8')+Raster(r'remap9')+Raster(r'remap10')+Raster(r'remap11')+Raster(r'remap12')+Raster(r'remap13')

now how do I execute this command and save it to new Raster file?

b=Raster(reclass1)
b.save("blah")

Traceback (most recent call last):
  File "R:\Data\LUCZ_2017\ConservationPriorityAreas\cons.py", line 36, in <module>
    Raster_Calculator(env2)
  File "R:\Data\LUCZ_2017\ConservationPriorityAreas\cons.py", line 34, in Raster_Calculator
    b=Raster(reclass1)
RuntimeError: ERROR 000732: Input Raster: Dataset Raster(r'remap1')+Raster(r'remap2')+Raster(r'remap3')+Raster(r'remap4')+Raster(r'remap5')+Raster(r'remap6')+Raster(r'remap7')+Raster(r'remap8')+Raster(r'remap9')+Raster(r'remap10')+Raster(r'remap11')+Raster(r'remap12')+Raster(r'remap13') does not exist or is not supported

update

reclass1 = [Raster(x) for x in arcpy.ListRasters() if x[:5]=='remap' if int(x[5:]) < 14 ]
print reclass1

print sum(reclass1)

use the regular python sum tool??

Traceback (most recent call last):
  File "C:\Users\Desktop\abid.py", line 19, in <module>
    print sum(reclass1)
RuntimeError: ERROR 999998: Unexpected Error.


texts: Copying attribute table to Excel in ArcGIS Pro?

texts:
In ArcMap I can select all rows in the attribute table, right click one of the grey boxes, copy and then paste into excel:

There is no Copy option in Pro.

How can this be done in ArcGIS Pro (without having to export the table)?


texts:
In ArcGIS, is it possible to change the polyline symbology to show vertices? Or do I need to extract the points via Feature Vertices To Points? This becomes a pain if I want to edit the source polygon/polyline features.

Here is what I'd like to see (this is from GMS), as I can clearly see start/end nodes and interior vertices.

I found out that pressing the V key while editing will show vertices. But it is only useful while editing and pressing V.


texts: Raster Subtraction using rasterio or gdal

texts:
I'm beginner to python. Each directory has two rasters named A.tiff and B.tiff. I need to subtract B.tiff from A.tiff and store the resultant file in the same directory. The directory paths are named by their date of creation(ex: /home/hash/images/20210121T123308/*A.tiff"). In a similar way I have directories created for each week from January with files A and B.
Below is my example code to fetch A and B's path
import glob
A_path = glob.glob("/home/hash/images/**/*A.tiff", recursive = True)
B_path= glob.glob(("/home/hash/images/**/*B.tiff", recursive = True)

Note: Both the file's metadata are same.


texts: How to export QGIS Contours to be compatible with 12D Software

texts:
I have processed my aerial drone survey images with Ground Control Points using AgiSoft. My accuracy and results were good.
I have exported the DEM into QGIS and generated the Contours with labels etc... I'm fairly certain the contours I have generated include a 3D vector with Z coordinates.
But now I need to know how to get this data into a format that my client (who used 12D Software and CAD) can use?
I guess the XYZ coordinates and contours are very important, but I am unsure how to export this in a way that is helpful for them.
Thanks!


texts: How to Unselect Features with QGIS tool "selectbyexpression' and the ">" Operator

texts:
First of all I've successfully used the "selectbyexpression" tool to UN-SELECT points by setting its 3rd parameter to "2" in the following script:

Input = QgsMapLayerRegistry.instance().mapLayersByName('Central_Region_civic_pts')[0]
clear_expression = '"STRNAME" = '+"'"+"Main"+"'"
processing.runalg('qgis:selectbyexpression', Input, clear_expression, 2)

But now I want to write an expression that will unselect EVERYTHING, not just a particular street by name (e.g. "MAIN").  This will be useful near the start of my script to wipe the slate clean, just in case anything is already selected, before my script goes on to make the selections I actually want.

Since my 'Input' file has a "PID" column ("Property ID", all numbers) I figured something like the following would work:

clear_expression = '"PID" > 0'

But it doesn't.

QUESTION

I suspect the problem is with the ">" symbol.  Can anyone confirm this, and suggest a workaround or another tool to do the job?

I also suspect there's nothing unique about the ">" symbol.  I'll bet there are more operator symbols than ">" that also don't work with this tool because the "Select By Expression" GUI at the top menu of QGIS only has 9 operator options at the top and ">" isn't one of them.  Neither is "<".  And there's no option for "does not equal" either, etc, etc...


texts: Enable Apply Edits With Global Ids on ArcGIS Feature Services

texts:
I have published a Feature Service on our ArcGIS Server 10.4.1 and in support section of the Feature Service URL I am seeing that the

  ApplyEdits With Global Ids: false

can you please let me know why this is false? and how can I update it to True?


texts: Converting bunch of TAB and TIF files into single GeoTIFF using GDAL?

texts: I have a number of TAB and TIF files. They form a single coverage. What I want to do is to convert all of them into a single GeoTIFF file, using GDAL tools, but I do not know how. I looked through dozens of threads, tutorials and etc., and it seems like GDAL supports TAB and TIF files, though I could not find any real example.


texts: Intersecting 2 layers without creating duplication using QGIS?

texts:
I have 2 polygons that i wish to combine into one layer, in QGIS.

on the first layer (neighborhoods) i have bigger areas (that divide a city to 7 major areas, or smaller neighborhoods).

that i wish to intersect with a second layer (statistical areas) that is built of smaller polygons- this layer holds population numbers and other data.

when I perform intersect, from the second layer i get duplication of the areas, so that each smaller area can show in one or two of the larger areas.
i want to perform this combination, while either:

1) setting a rule that each of the smaller area will be associated only with one of the bigger areas, the one that the majority of the polygon is in

or

2) create the same intersection with the duplication but that an attribute of my choice (population number) will be divided relatively to the size of the area that is relevant to the bigger area (under the assumption that the population is equally distributed)


texts: Extracting vector/raster data for CSV data points

texts:
I am using QGIS and I have a CSV file of lat and long coordinates which represents presence points of a species in the USA. I also have vector (the state extents/names) as well as raster information (NDVI, bioclim, elevation, etc.). I would like to be able to extract data from both the vector and rasters which is associated with each coordinate. For example, so I can say X amount of coordinates were in Kentucky, and X points were at an elevation greater than 500 m.

At this point I would just like to ask if anybody can please suggest any QGIS tools or plugins that are able to do this?