texts: Replacing multiple polygon features by their shared center point in GRASS GIS

texts:
I have three polygon features of the same category in a larger vector layer from which I would like to create a center point (mean coordinate) that holds the same attributes as the three polygons.
I would then like to delete the original polygon features and replace them with the new point in the larger vector layer.

Running the below commands finds me the center point, which, following the answer by Micha, can be connected to the attribute table of the polygons to hold the same data. However, during these steps the character columns are redefined as text columns, making it so that it can't be patched with the original vector layer as 'Length of string columns differ'. I thus need a way to avoid this transformation, or a way to transform the column-types back again. I will look into this last point, which I assume can be done by defining a new table and copying over the data - although that sounds like it should have a better solution.

# First extract the three polygons that have cat 1
v.extract input=vector_layer cats=1 output=polygons

# Content of the attribute table
v.db.select polygons
cat|precision|method_con|measure_me|survey_yea|event|elev_minimum
1|3|tpDmg|gnss|2017||101.230186462402

# With corresponding column names and types
v.info -c polygons
Displaying column types/names for database connection of layer <1>:
INTEGER|cat
INTEGER|precision
CHARACTER|method_con
CHARACTER|measure_me
INTEGER|survey_yea
CHARACTER|event
DOUBLE PRECISION|elev_minimum

# Then find the center points of the individual polygons
v.centerpoint input=polygons output=poly_center

# These automatically hold an identical attribute table as the polygons
v.info -c poly_center
Displaying column types/names for database connection of layer <1>:
INTEGER|cat
INTEGER|precision
CHARACTER|method_con
CHARACTER|measure_me
INTEGER|survey_yea
CHARACTER|event
DOUBLE PRECISION|elev_minimum

v.db.select poly_center
cat|precision|method_con|measure_me|survey_yea|event|elev_minimum
1|3|tpDmg|gnss|2017||101.230186462402

# Then find the center point of the three polygon center points
v.centerpoint input=poly_center output=centerpoint

# This point does not have an attribute table.
# Add attribute table and join with that of the poly_center points
v.db.addtable centerpoint
v.db.join centerpoint column=cat other_column=cat other_table=poly_center

# The data itself is the same
v.db.select centerpoint
cat|precision|method_con|measure_me|survey_yea|event|elev_minimum
1|3|tpDmg|gnss|2017||101.230186462402

# However, the 'character' columns (all of length 254)
# have been redefined as 'text' columns of length 1000.
v.info -c centerpoint
Displaying column types/names for database connection of layer <1>:
INTEGER|cat
INTEGER|precision
TEXT|method_con
TEXT|measure_me
INTEGER|survey_yea
TEXT|event
DOUBLE PRECISION|elev_minimum

# Delete the three polygons from the larger vector layer
db.execute sql="DELETE FROM vector_layer WHERE cat = 1"

# Attempt to patch the centerpoint to the larger vector layer
v.patch -e input=vector_layer,centerpoint output=new_vect
ERROR:Length of string columns differ

I am using GRASS GIS 7.6.1 on Windows 10. I am also open to a complete reworking of the steps taken, as they certainly might not be the best.


texts: QGIS CRS Coordinate Reference System -- Generated vs. Specified -- error levels?

texts:
Below image appears to indicate that the vector layer data was identified as
NAD_1983_HARN_StatePlane_Florida_North_FIPS_0903_Feet, which is a well known (to QGIS) CRS. However when the layer was added it was automatically converted to WGS 84. (I presume WGS is an excellent "universal" CRS?)

Is the conversion algorithm inside QGIS to WGS just (about) as accurate to telling QGIS the specific source CRS (e.g. NAD_1983_HARN_StatePlane_Florida_North_FIPS_0903_Feet) ? Is there great advantage to specifying the CRS vs. just allowing the automatic conversion to WGS to occur?

I will likely add a OSM layer (not sure what CRS that will be) and may need to add other layers from the U.S. Census, Google Maps, etc. Wondering if I should pay close attention to the CRS specifics of each layer, or just rely on QGIS to convert to WGS and call it a day?

I understand in general the basics in 7.1 Lesson: Reprojecting and Transforming Data, but remain uncertain on the specifics.


texts: OSRM - trucks and cars routes should differ

texts:
I am working on logistic software that includes building routes by points given. So I've got my local OSRM server and it's working fine.

But now I faced the problem of dividing vehicles into cars and trucks. Some areas are conditionally closed for the truck passage, so I need the truck route to be built bypassing that areas.
So I'm wondering how do I get this?

I must admit I'm new to OSRM and OSM data, so I'm thinking about following steps:

If forbidden area is not in OSM data, should I add it to my own data (like in this question)?
If an area exists or I've added it, can I set up OSRM profile to consider the area as "inaccessible for the trucks"? Should I re-extract map data after that?
Can I pass the param "has acess to forbidden areas" when quering for a route?

Am I on the right way and is my goal reachable?


texts:
I have a rasterA (whole Scandinavia, with values from 0-255) and a rasterB (in blue, a small part of south of Norway, entirely inside rasterA). I want to query the pixels of rasterA which are inside rasterB.

Will this query return only those pixels which overlap, or will it return values for the whole of rasterA, since rasterA and rasterB overlap?

Query:
SELECT A.rid, pvc FROM rasterA AS A, rasterB AS B, ST_ValueCount(A.rast,1) AS pvc WHERE st_intersects(A.rast, B.rast);

rasterA(white) and raster B (blue):

Also, I'm importing the raster into PostGIS like so:

raster2pgsql \
    -s 3857 -I -C -Y \
    -t ${TILESIZE} \
    -l 4,32,128 \
    ${RASTERFILE} $TABLENAME |
    psql -q --set ON_ERROR_STOP=1

Can this have anything do to with the faulty results I'm getting?


texts: Placeholder Text in search box - Cartodb

texts:
I am trying to add placeholder text to the Cartodb searchbox using the following DOM code

document.getElementsByClassName('text')[0].placeholder='Search Your Constituency';

It actually does not seem to work.


texts:
I am making SVG files with online servers like "online-convert" or picsvg and with the print composer.

Me problem is that QGIS is not working well with files that are not located in the SVG folder. This folder is a sub folder in the QGIS downloads library and it is taking relatively long time to open this folder in order of saving new SVG files.
Is there a way for an easy and fast way to locate this folder?


texts: Selecting multiple date ranges using Select By Attribute?

texts:
Using ArcGIS Desktop 10.6.1. I'm interested in selecting multiple date ranges using Select By Attribute.

I've verified that each date range works. For example Date >= '2015-04-15' AND Date <= '2016-04-15' works great.

But I would like to also select that same date range for 2016-2017. I get an error if I continue with AND for the next range.


texts: Why do points shift when reprojected from WGS84 to a state plane coordinate?

texts:
We have a Trimble GPS unit we use to collect system assets.

We collect the data in WGS84 and differentially correct it, and then reproject it into our State Plane Coordinate System.

The resulting points always shift the same distance in the same direction after the reprojection.

Any idea what could cause this?


texts:
I am looking at bottom friction (Manning's n) values in relation to storm surge as large volumes of water flow over a surface.

In order to calculate total head loss occurring as water passes from the original water's edge to the final flooding extent, I must somehow figure out a way to account for cumulative losses. Is there a way I can create a raster map containing the cumulative (or average) Manning's n values over which the water has traveled?

In other words, with the accumulation of friction starting from the original water's edge, the first cell would contain its own Manning's n value (I will call this n1). Moving away from the water's edge, the next cell will be n1+n2. The next cell away from the water will be n1+n2+n3 and so on.

Below is a script I have created for making a cumulative Manning's n map for a small (121-cell) sample raster of Manning's values. The sample raster is named manningWithWater and this is done using r.buffer. The buffer starts from the raster called falseWater, which is artificial water values I have added to the Manning's raster. The script works but is inefficient; this took 22 seconds on my tiny sample raster and would use up quite a bit of memory while processing a regularly-sized raster.

import grass.script as grass

def main():
    x = range(100,700,50)
    numString = ""
    for num in x:
        numString = numString + str(num) + ","
    numString = numString[:-1]
    grass.run_command('r.buffer',
        overwrite=True,
        input="falseWater@PERMANENT",
        output="falseWaterBuffer",
        distances=numString,
        units="feet")
    totalBuffers=grass.read_command('r.describe',
        map="falseWaterBuffer")
    maxBuffer=int(totalBuffers.split("-")[-1])
    finalExpression=""
    for n in range(2,maxBuffer+1):
        grass.mapcalc(
            "$output=if($buffer==$bandNumber,$manningWithWater,null())",
            overwrite=True,
            buffer="falseWaterBuffer@PERMANENT",
            output="band"+str(n),
            bandNumber=n,
            manningWithWater="manningWithWater@PERMANENT")
        grass.run_command('r.grow.distance',
            overwrite=True,
            input="band{0}@PERMANENT".format(n),
            value="band{0}Grown".format(n))
        grass.mapcalc(
            "$output=if($buffer>=$bandNumber,$bandGrown,0)",
            overwrite=True,
            output="band{0}Add".format(n),
            buffer="falseWaterBuffer@PERMANENT",
            bandNumber=n,
            bandGrown="band{0}Grown@PERMANENT".format(n))
        finalExpression=finalExpression+"band{0}Add@PERMANENT+".format(n)
    finalExpression=finalExpression[:-1]
    grass.mapcalc(
        "cumulativeN={0}".format(finalExpression),
        overwrite=True)

if __name__ == '__main__':
    main()


texts: How to obtain the difference between two time fields in hh:mm:ss format?

texts:
I have to 2 date fields in a table, ALARMTIME and ARRIVETIME.  I want to acquire the difference between the two in a new RESPONSE_TIME field in ArcGIS.
The only successful field calculation I could get is the answer in seconds using:
DateDiff ("s", [ALARMTIME], [ARRIVETIME] )
But if I have 8:05:00am - 8:02:32am, I'd like to see 00:02:28 or 02:28 or 2:28 instead of 148 seconds.  Please help.  My response time field is a double.  I couldn't get any calcs to work using a date field.


texts: Supervised classification: Why is the label built-up area not showing in the graph but instead the number 8?

texts:
Once exporting my work the graph seems to show the number 8 instead of the label built up area. What is the problem?
This is my code:
var image = ee.Image(ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
    // Filter to get only images under the region of interest.
    .filterBounds(roi)
    // Filter to get only one year of images.
    .filterDate('1986-01-01', '1986-12-30')
    // Sort by scene cloudiness, ascending.
    .sort('CLOUD_COVER')
    // Get the first (least cloudy) scene.
    .first());

Map.addLayer(image, {bands: ['B5', 'B4', 'B3'], max: 0.5, gamma: 2}, "L8 Image");
// Merge points together
var newfc = Water.merge(naturalwoodlands).merge(plantedforest).merge(builtuparea).merge(cultivatedland).merge(grassland).merge(barrenland).merge(saltmarsh);
print(newfc, 'newfc');
// Select the bands to be used in training
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'];

// Sample the input imagery to get a FeatureCollection of training data.
var training = image.select(bands).sampleRegions({
  collection: newfc,
  properties: ['landuse'],
  scale: 30  // should reflect the scale of your imagery
});
// Make a classifier and train it.
var classifier = ee.Classifier.randomForest(10).train({
  features: training,
  classProperty: 'landuse',
  inputProperties: bands
});
// Classify the input imagery.
var classified = image.select(bands).classify(classifier).aside(print);

// Define a palette for the Land Use classification.
var palette = [
'#0a9917', // naturalwoodlands(0)
'#00f7f7', // water(1)
'#fffcb0', // cultivatedland(2)
'#d7d8ff', // builtuparea(3)
'#1aff0b', // plantedforest(4)
'#e766ff' , // grassland(5)
'#1612f7', // barrenland(6)
'#a054ff',// saltmarsh(7)
];
// Display th#ff4b38e classification result and the input image.
Map.setCenter(23.054, -34.0482);
Map.addLayer(classified, {min: 0, max: 6, palette: palette}, 'Land Cover Classification');
// Get a confusion matrix representing resubstitution accuracy.
var trainAccuracy = classifier.confusionMatrix();
print('Resubstitution error matrix: ', trainAccuracy);
print('Training overall accuracy: ', trainAccuracy.accuracy());

// var ghnp = ee.FeatureCollection('ft:1Kh5uxAdbnuF77vydlLPCAOQ-A7HgLxyYFyI4QnnT', 'geometry');
// Map.addLayer(ghnp, null, 'roi');

var areaChart = ui.Chart.image.byClass({
  image: ee.Image.pixelArea().divide(1000000).addBands(classified),
  classBand: 'classification',
  region: geometry,
  scale: 30,
  reducer: ee.Reducer.sum(),
  classLabels: ['water', 'naturalwoodlands','plantedforest','builtuparea','cultivatedland','grassland','barrenland','saltmarsh']
}).setOptions({
  hAxis: {title: 'Classes'},
  vAxis: {title: 'Area sqkm'},
  title: 'A histogram showing area by class',
  series: {
    0: { color: 'blue' },
    1: { color: 'green' },
    2: { color: 'brown' },
    3: { color: 'gray' },
    4: { color: 'yellow' },
     5: { color: 'purple' },
      6: { color: 'pink' },
       7: { color: 'black' },
  }
});
print(areaChart);
var exportAccuracy = ee.Feature(null, {matrix: trainAccuracy.array()})

// Export the FeatureCollection.
Export.table.toDrive({
  collection: ee.FeatureCollection(exportAccuracy),
  description: 'exportAccuracy',
  fileFormat: 'CSV'
});
// Get a confusion matrix representing resubstitution accuracy.
print('RF error matrix: ', classifier.confusionMatrix());
print('RF accuracy: ', classifier.confusionMatrix().accuracy());
Export.image.toDrive({
 image: classified,
   description: 'LC_2018',
   scale: 30,
   region: geometry,
maxPixels: 130000000000,
 });
var exportAccuracy = ee.Feature(null, {matrix: trainAccuracy.array()})

// Export the FeatureCollection.
Export.table.toDrive({
  collection: ee.FeatureCollection(exportAccuracy),
  description: 'exportAccuracy',
  fileFormat: 'CSV'
});


texts: What is the best plain text format for georeferenced raster data?

texts:
I have to translate satellite images in a file in .txt format or another, such as Excel or .csv, in which there are pixel values corresponding to (longitude,latitude) couples. What is the best file format to store this (in terms of dimension and limits of the file because pixels number is in the order of 10^5)?

In your experience, is a matrix format representation better or a columns representation (one column is for longitude, one for latitude and one for value, so every row of the text corresponds to a pixel)?


texts: I want to load multiple raster data files into a QGIS project. In total the are about 2,5 GB in size. Is there a way to pre-unselect the layer prior to the loading process to not crash the whole system?


texts: arcpy.geometry __geo_interface__ and AsShape() function: loss of precision and holes

texts:
I'm serializing my arcpy geometries as geojson so that I can 'hydrate' them back as geometries later and I'm having 2 problems in the cycle.:

PROBLEM 1: Precision

    R0 = arcpy.SearchCursor(self.shpTest, "FID=0").next().getValue("Shape")
    geojson = R0.__geo_interface__
    R1 = arcpy.AsShape(geojson)
    self.assertTrue(R0.equals(R1)) <<< THIS FAILS

If I check the string representation, the coordinates have slightly changed:

    geojson2 = R1.__geo_interface__
    print geojson
    print geojson2

    {'type': 'Polygon', 'coordinates': [[(442343.5516410945, 4814166.6184399202), (442772.17749834526, 4811610.7383281607), (441565.67508534156, 4811499.6131059099), (440772.50052100699, 4814184.7808806188), (442343.5516410945, 4814166.6184399202)]]}
    {'type': 'Polygon', 'coordinates': [[(442343.55169677734, 4814166.6185302734), (442772.17749023438, 4811610.73828125), (441565.67510986328, 4811499.6130981445), (440772.50048828125, 4814184.7808837891), (442343.55169677734, 4814166.6185302734)]]}

PROBLEM 2: Holes
If the polygon has holes, geo_interface generates an error:

    R0_WithHoles = arcpy.SearchCursor(self.shpTest, "FID=0").next().getValue("Shape")
    geojson = R0.__geo_interface__  <<< generates this ERROR:

    File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\geometries.py", line 68, in __geo_interface__
        return {'type': 'Polygon', 'coordinates': [[(pt.X, pt.Y) for pt in part] for part in self]}
    AttributeError: 'NoneType' object has no attribute 'X'

Any ideas on how to solve these problems?


texts: r.terraflow fails on Mac (rusage/gettimeofday: Bad address)

texts:
I am trying to get flow lines, flow accumulation and stream segments form a raster DEM. The DEM is about 2000x3000 cells.

I loaded this into GRASS (v7.3, running on Mac OS El Capitan) using a .prj projection definition because it is in a non-EPSG CRS. I set the region to fit the raster and match the cell alignment.

When I run r.terraflow I get a warning:

WARNING: raster Rush_Creek_DEM is of type DCELL_TYPE --precision may be lost.

The command ends with:

rusage/gettimeofday: Bad address

Any ideas how I can get this to work?

The entire output using the --verbose flag, shown below:

r.terraflow --verbose elevation=Rush_Creek_DEM@PERMANENT filled=Rush_Creek_DEM_fill direction=Rush_Creek_DEM_flowdir swatershed=Rush_Creek_DEM_watershed accumulation=Rush_Creek_DEM_flowaccum tci=Rush_Creek_DEM_tci
Memory manager registering memory in MM_WARN_ON_MEMORY_EXCEEDED mode.
cell Rush_Creek_DEM header compatible with region header
Elevation stored as FLOAT (4B)
WARNING: raster Rush_Creek_DEM is of type DCELL_TYPE --precision may be lost.
Region size is 2129 x 3212
STREAM temporary files in /Users/brendan/grassdata/Rush_Cree
k/PERMANENT/.tmp/MBP.local/13929.0  (THESE INTERMEDIATE
STREAMS WILL NOT BE DELETED IN CASE OF ABNORMAL TERMINATION
OF THE PROGRAM. TO SAVE SPACE PLEASE DELETE THESE FILES
MANUALLY!)
Reading data from <Rush_Creek_DEM> to stream </Users/brendan/grassdata/Rush_Creek/PERMANENT/.tmp/MBP.local/13929.0/STREAM_baaaaX>
total elements=6838348, nodata elements=145083
largest temporary files:
FILL: 365.21M (382947488) [6838348 elements, 56B each]
FLOW: 510.66M (535461200) [6693265 elements, 80B each]
Will need at least 1021.31M (1070922400) space available in /Users/brendan/grassdata/Rush_Creek/PERMANENT/.tmp/MBP.local/13929.0
EMPQUEUEADAPTIVE: starting in-memory pqueue
EMPQUEUEADAPTIVE: available memory: 297.99MB
EMPQUEUEADAPTIVE: desired memory: 297.99MB
sz_stream: 263232 buf_arity: 200 mm_overhead: 8437120 mm_avail: 312465395.
EMPQUEUEADAPTIVE: memory overhead set to 8.04626MB
EMPQUEUEADAPTIVE: pqsize set to 38003534
rusage/gettimeofday: Bad address


texts: CartoDB-Picture not showed inside Infowindow

texts:
I am trying to customize my infowindow in cartodb but the picture is not showed.
I have the same picture uplooad into google drive and dropbox in two differents formats *.jpeg and *.png but it does not work....

Any idea why? Maybe I must upload in another public place....


texts: Alternative to 'for' loop using list of strings in Google Earth Engine

texts:
I am new to Google Earth Engine code.
Currently, I have code that runs various functions (NVDI, greenness index from Landsat data) for a single state, but for several years (using a 'for' loop for the years) - it generates a NDVI and other factors for 1997, 1998,...2019.  I run one state, change the name, and run it again.  That's cumbersome, and since I have many states to do, so I'm hoping to automate it.
example of current code:
var startYr = 1997;
var endYr = 2019;
var roiName = 'TN'
var states = ee.FeatureCollection('my_assets/myshapefile'); // shapefile with states
var roi = ee.Feature(states.filter(ee.Filter.eq('STUSPS',roiName)).first()); // State for ROI

for (var yr = startYr; yr <= endYr; yr=yr+1){
      blah blah blah function to calculate Landsat band metrics for the roi for each year blah blah blah
}

Now I want to run this but not have to change the code for each state I want to run.
So I have a list of roi's:
Example:
var states = ee.FeatureCollection('my_assets/cb_2018_CONUS_500k'); shapefile with all U.S. states
 var filtered = ee.Filter.inList('STUSPS', ['AR','MS','LA','MO','TN']); //list of states to process

I guess I'm stuck on how to "properly" do the 'for' loop alternative for the list of roi's.
I've snipped a bunch of the code out so hopefully this makes enough sense.


texts: GIS-based reporting and maintenance tool

texts:
Background

Conceptually, I am looking for a tool, or suite of software tools that can be readily integrated to help browse data using a GIS-based map (such as Google Maps). The proposed architecture resembles, but is not limited to, the following:

The technologies shown above are not fixed. For example, if OpenGEO can eliminate many of the components while adding others that are not listed (such as authentication and authorization), that would be most helpful.

Essentially, the architecture boils down to:

Map -> Filter sites -> Select sites -> CRUD

And:

Map -> Filter sites -> Select sites -> Summary Reports

The following sections describes these two problems in more detail, using a hypothetical (broad-brushstroke) example.

Problem - CRUD

In the first screen, the user is prompted to filter by some arbitrary criteria (team names, regions, skills, equipment, and more) and select some sites:

The sites that match the filtered options are shown on the map and the user is then presented with a number of options that can be used to edit data:

If the user chooses Personnel, the system provides a CRUD screen:

The CRUD screens can be fairly complex, having simple code tables to drive them, such as:

Problem - Report

The second part of the problem is very similar to the first, but instead of visiting CRUD screens, the user can choose to generate reports:

Note that the user has already selected the sites for reporting. If the user selects the equipment option, a report is generated with summary data:

Question

What off-the-shelf software tool(s) can provide most (if not all) of these features?

Note: The interaction need not be exactly as shown, but the main use cases must be possible.

Use Case 1 - Populate Database

Define regions on a map.
Put sites [building locations] on a map.
Insert data (for each site).

Use Case 2 - Selection

System shows a map.
User filters by region and other criteria.
User selects sites.
System provides options (maintenance or reports).

Use Case 3 - Maintenance

User selects maintenance.
System presents complex CRUD.
User can add, edit, or delete data.

Use Case 4 - Reports

User selects report.
User selects type of report to generate.
System generates summary report.


texts: How to add a GeoServer WMS using OpenLayers-3?

texts:
I am unable to figure out why I can't add WMS in OL3? Tiger is a default layer in Geoserver.

Below my code.

    Open Street Map

&nbsp;


texts: Is it possible to take advantage of lowercase index on column using CQL filter on WFS?

texts: I am querying a WFS housed in geoserver - the data itself lives in a postgresql database. I would like to run case insensitive queries against the WFS and see there is an option to specify strToLowerCase on the field I am querying, but it doesn't appear to take advantage of the lowercase index I have created in postgresql. The table I am querying has several million records, and when using the strToLowerCase I run into memory issues. Is there a way I can do a case insensitive query on the data in the WFS and take advantage of the lowercase index that exists?


texts: add definition query to layers in a group

texts:
I'm trying to develop a python code that will apply the same definition query to all layers in a specified group.  I currently have a working code that will apply a definition query to all layers in the Table of Contents, but unfortunately I need it to apply the syntax to a specific group of layers. Any ideas?

I've seen this post:

Create a definition query for groups of layers: Is it possible?

But when I set it up specific to my mxd project, nothing happens.

import arcpy

#Variables to form query syntax
#-------------------------------------------
#field in Attribute table
field = 'Country \n'
queryField = '"%s"' % field.strip()   #add double quotes field (for query syntax)
#value in specified field
value = "'Canada'"
#concatenate query syntax
queryStr = str(queryField) + "=" + str(value)
#--------------------------------------------
mxd = arcpy.mapping.MapDocument("CURRENT")

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DEFINITIONQUERY"):
        lyr.definitionQuery = queryStr
arcpy.RefreshActiveView()
del mxd


texts: Georeferenced Tiffs displaying in the wrong location & scale when updating ArcGIS 10

texts:
I'm working on a project which involves georeferencing a large number of .tiff files to digitize data from them. All the tiffs which have already be georeferenced where done with ArcMap 10.0 without any service packs installed. In an effort to solve an issue with the pixels inverting and a certain zoom level, I've been experimenting with upgrading service packs. Each service pack I install causes some, but not all, of the previously georeferenced tiff files to display at a significantly smaller size, and in the incorrect location. The location change may just be do the the scaling problem though, since it is still near the correct location, I'm not positive though. Also, when I open the same map document on the other computer in the office running ArcMap 10.0 without any service packs, the tiffs are located where they should be and the proper size. I'm currently running ArcMap 10.0 with Service Pack 5 installed, and I've also installed the Raster Format Patch, neither of which have helped. I've done a few searches, and have not been able to find any solutions for this issue. Has anyone else run into this problem?

Thanks,
Chris


texts: concurrent use of in_memory workspace in arcpy

texts:
i am using a python script that's eliminating small polygons from a feature class. this script creates temporary feature classes with fixed (hardcoded) names in the "in_memory" workspace. example:

arcpy.CopyFeatures_management("in_fc","in_memory/item_a")

it's a fairly reusable script, so i use it at many points in a geoprocessing chain.

am i in danger, that the same script called twice from different processes (at the same time) will have a conflict at, for example "in_memory/item_a" ? or are the "in_memory" workspaces seperate for each process?

i am using ArcInfo 10.0.


texts: I currently have a 5m DEM with 6000 rows and columns. I need to keep the spatial resolution but reduce the number of rows and columns to 1000. I have used the basin tool to separate the DEM and extract by mask to split up the DEM into smaller sections the size of the basins but the rows and columns are still 6000. Does anyone know anyway round this? I think I may be able to use Environments to do this but unsure. I am using ArcMap.


texts: Run OTB commands in OSGeo4W Shell scripts

texts:
I am trying to run shell scripts with Orfeo-toolbox commands, and I'm new to all this.... (I'm working in Windows 7).

I can type

otbcli_ComputeImagesStatistics -il Imagery/P1_10cm1.tif Imagery/P2_10c1.tif -out ImageAnalysis/Stats12_JUNK3.xml

directly into the OSGeo4W shell window and things work fine. Now I'd like to put more commands in a script and run them all at once.

I'm following the advice from this question.

If I save that text in a file called TryShell.sh
and enter in shell window:

sh TryShell.sh

I get the error TryShell.sh: line 1: otbcli_ComputeImagesStatistics: command not found

Why does the otbcli command work in the shell but not in the shell script?


texts: How to add value to field according to other fields by using where condition in QGIS3.0?

texts:
I have a dictionary created from python. I want to update my attribute table accroding to that dictionary in QGIS3.0:

dict = {1:4, 2:5, 3:8, 4:9, 5:12, 6:17}

I want to see that values according to something like this condition:

UPDATE su_kutlesi = 1 WHERE Subbasin = 4
UPDATE su_kutlesi = 2 WHERE Subbasin = 5
UPDATE su_kutlesi = 3 WHERE Subbasin = 8
UPDATE su_kutlesi = 4 WHERE Subbasin = 9
UPDATE su_kutlesi = 5 WHERE Subbasin = 12
UPDATE su_kutlesi = 6 WHERE Subbasin = 17

I looked CASE from conditionals in qgis:fieldcalculator. It seems okey for me but I can not use that. Also, I do not want to insert loop my code but if neccessary, it is okey.Here is my attribute table:


texts: How to use a common GeoWebCache for multiple Geoservers?

texts: I had multiple Geoserver's (5) in different systems one is back up for another with same port and same layers are published . I know that every Geoserver is integrated with Geowebcache. Now I want to turn off individual GeoWebCache on every system and to install a new GeoWebCache  for all geoserver's on any system.How to do?


texts: Check if raster (part of) lies within a pixel - python/gdal

texts:
I have a list of the corner-coordinates of different grid (c1-c4). First i sorted these coordinates to ensure that c1 is the upper left, c2 is the upper right, c3 is the lower left and c4 is the lower right coordinate. My second input is an empty grid with a given cellsize.

Now i want to get through every pixel in the empty grid and check with the four coordinates if they within the pixel of the empty grid.

It would be nice if somebody could help my. First i want to know how the get the x- and y-coordinate of every pixel in the grid.
Once the grid (blue) is inclined also the edge regions (grey) are excluded (see picture).


texts: Creating curved annotations for polygons and lines in QGIS

texts:
In general, my problem is that: I can`t display annotations using the option "Using the perimeter (curved)" in QGIS. I tried to solve it in two ways:

At first was an attempt to use just a simple placement option "Using the perimeter (curved)" in the polygon layer, but when I chose, I lost annotation on the map.
Next step was when I tried to convert polylines into a line and then display my curved annotation, it also didn't work, as my result wasn't a line drawn through the centre of the polygon, but a frame around the polygon?


texts: "The tool not licensed" ERROR on a geoprocessing service on ArcGIS Server 10.1

texts:
I have a python script which runs well as geoprocesing tool. When published on ArcGIS Server 10.1 as geoprocessing service, however, it throws the following error:

Unable to complete operation.
Error executing tool.: Failed to execute. Parameters are not valid. ERROR 000824: The tool is not licensed. Failed to execute (Kriging). Failed to execute (NewCuringModelExpPy). Failed to execute (New_Curing_Model_ExpPy).

The script runs without any error on ArcCatalog with Spatial Analyst extension authorized. The code also has got proper license checked out as below:

arcpy.CheckOutExtension("spatial")

I don't think ArcGIS Server requires Spatial Analyst extension authorized, does it? Otherwise, is it to do with any ArcGIS Server configuration?

By the way, I am running on ArcGIS Server 10.1 Advanced so Spatial Analyst extension should have been part of it.


texts: How to find initial values for fit.variogram

texts:
I wish to estimate a variogram for spatially distributed price data in Tanzania. I'm new to both spatial statistics and gstat, and have a couple questions.

(1) I can SEE spatial patterns in price, and moran's I is about 0.2, so I know there is spatial correlation. But when I plot my variogram using the code below, semivariance is flat across distance. (Screenshot in dropbox link below, as well as the tzprice1 shapefile for replication.) Why would this be? Perhaps I need to specify smaller intervals, and a smaller range? And if so, how?

tzprice1_v <- tzprice1[!is.na(tzprice1@data$prodp), ]
pr.v <- variogram(prodp ~1, tzprice1_v)
plot(pr.v)

(2) I am unable to estimate the variogram using fit.variogram, and I'm guessing that this is because the underlying data don't show the right shape. Is this correct? But just in case, I'm using the code below... is it correct to simply guess the sill, nugget, and range starting values in this way, from the plot that I did above?

fit.variogram(pr.v, vgm(psill=.125,"Exp",range=200,nugget=.1))

Shapefile and screenshot here:
https://www.dropbox.com/sh/oad1xfcpcugkulh/AACec_FuUZqCHR9-Jicohgm9a?dl=0


texts: Calculate distances from points to multiple country borders in QGIS

texts:
I have a shapefile containing several thousand GPS locations located in multiple countries and another shapefile with country polygons (or alternatively also one with country border polylines).
How can I calculate the shortest distance from each GPS location (point) to the borders of all countries adjacent to the country the point is located in?
Is there even a way to achieve this in QGIS?
Example of what I would like to calculate for every point:

I managed to get the distance to the (single) nearest border using the NNJoin plugin, but I have trouble identifying all adjacent borders and calculating all distances.
Ideally, the data structure would look like this, assuming that the distance is only reported for adjacent countries:

If this is difficult to achieve for adjacent countries only, a simple distance matrix with the distances to all countries in the sample would also be fine, but I have not found a way to create one for polylines rather than points.


texts: Does QGIS 3.4 locks the row while editing the postgis data?

texts: I want to know if QGIS locks the row of a postgis table for editing  or can this behavior be force implemented? Suppose there is another user who wants to edit the same row-- this must not be allowed.


texts:
I've installed mapnik3, python bindings, and postgis (and imported a pbf file into postgis), and can't seem to even make it generate a simple map as a png file.

Specifically - where can I find examples of a working style.xml file which will actually work with freshly imported OSM data. All tutorials that I can find are related to mapnik 2.x.


texts: NDVI mask Error "image.normalizedDifference is not a function"

texts:
I'm new in Google Earth Engine and want to mask NVDI but get an error. How can I solve the problem?
This is my code: Link


texts:
I'm looking for the algorithms (the math) used by QGIS but I can't find any documentation about it

The ones I'm looking for are:

Minimum bounding geometry (convex hull)
Clip a polygon by a polygon

Does anyone knows where I can find them?


texts: ESRI JavaScript API: Feature layer from x,y events

texts:
I am trying to create a Feature Layer on the ESRI JavaScript client using the Map Service Feature layer that was created on the map document as an XY event theme (ArcMap -> File -> Add Data -> Add XY Data).  The data source for this is a 'view' on the SQL server and therefore this view doesn't have an Object ID!
The feature Layer thus generated doesn't seem to load all the features!  Only the first fetaure is diaplayed although there are total of 7 features.


texts: Creating DEM from LAS file without using LAStools?

texts:
I want to create a ground model from a *.las file. The LAStools are available in QGIS but limited in use (generation of black diagonal in result file), so I can't use this.

In the GRASS tools I can find the tools v.lidar.edgedetection, v.lidar.growing and v.lidar.correction, to classify the las data. Unfortunately, v.lidar.edgedetection gives me a python error. So this isn't working either.

How could I proceed?


texts: SVM object based classification QGIS with Orfeo Toolbox. Model outputs incorrect/ invalid

texts:
I am trying to run an SVM object based classification on QGIS with Orfeo Toolbox. I have 5x Spectral Tif images of bands; blue,Green, red, Infra-red and Rededge, these were generated from a UAV survey using a multi-spectral camera.

My 'model' outputs from the training classification through 'Train images classifier' delivers multiple files such as shp. prj. dbf. etc.. but no singular model file. I have been thus far unsuccessful in running the classification. I have tried using all 5 spectral bands, then also just NDVI as well. Does anyone know how to save just one model file? or know where I may be going wrong here?


texts:
I have a .shp with points which have unique IDs based on which I want to connect them.
When I run the code, it creates empty line layers without geometry.

layer = QgsProject.instance().mapLayersByName("shap")[0]
features = layer.getFeatures()
hodList=[]
for feat in features:
    ids=feat['shape_id']
    hodList.append(ids)
unikatniId=list(set(hodList))

for i in unikatniId:
    expression = ('"shape_id" LIKE \'%s%s\''%(i,'%'))
    request = QgsFeatureRequest().setFilterExpression(expression)
    PointList = []
    v_layer = QgsVectorLayer("LineString?crs=epsg:4326", "line", "memory")
    pr = v_layer.dataProvider()
    line=QgsFeature()
    for f in layer.getFeatures(request):
        termino = QgsPoint(feat['shape_pt_lon'],feat['shape_pt_lat'])
        PointList.append(termino)
        line.setGeometry(QgsGeometry.fromPolyline(PointList))
        pr.addFeatures( [ line ] )
        v_layer.updateExtents()
        QgsProject.instance().addMapLayers([v_layer])


texts:
I want to use a free or at least a cheap alternative to ArcGIS Server that I can develop for in .NET.
MapGuide used to be such.

Is it still maintained?


texts: Reading BLOB data from an Esri GDB with OGR

texts:
Is it possible to read a binary data field (BLOB) from an Esri geodatabase using OGR? I'm using the OpenFileGDB driver and can extract the first few bytes, but then it truncates the data with "...".

>>> import ogr
>>> src = ogr.Open("data.gdb")
>>> src.GetDriver().name
'OpenFileGDB'
>>> layer = src.GetLayerByName("Photo__ATTACH")
>>> feature = layer.GetFeature(1)
>>> data = feature.GetField("DATA")
>>> len(data)
75
>>> data
'FFD8FFE000104A46494600010100004800480000FFE100584578696600004D4D002A0000...'

I'm not sure if this is a problem with the driver, or if I just need to pass it a particular flag?

FFD8 is the start of a JPEG file in HEX.

The equivalent in ArcGIS-speak is here: http://support.esri.com/technical-article/000011912


texts:
I am writing on an update of the qgis2leaf plugin.
It would be great to store the settings of the plugin UI (like which checkboxes are checked, which strings are used) in a separate file (like xml or sth. else). secondly it would be great to load such a file and set the UI with the parameters in this file.

Does anyone know a plugin for qgis that might already have such a functionality so I can use it as a template?


texts: Legend is not displaying all active features/layers in ArcGIS Pro using ArcPy?

texts:
I created a script tool that positions the north arrow, legend, scale bar, and title in a layout; and then zooms to a bookmark and exports the layout as a PDF. But when I go and look at the PDF, the legend does not display the graduated color symbology from the results of a summarize with-in, but will display one of my point layers. However, when I go and look at the layout in ArcGIS Pro the graduated color symbology is displayed in the layout legend.

Here is the custom script tool I'm using:

import arcpy

#input layer
lyr = arcpy.GetParameterAsText(0)

# input name of layout
p = arcpy.mp.ArcGISProject("CURRENT")
lyt = p.listLayouts("Layout_King")[0]

# Reposition the scale bar
scaleBar = lyt.listElements("MAPSURROUND_ELEMENT", "Scale Bar")[0]
mf = scaleBar.mapFrame
scaleBar.elementPositionX = mf.elementPositionX + 0.0
scaleBar.elementPositionY = mf.elementPositionY - 0.5

# Reposition the north arrow
northArrow = lyt.listElements("MAPSURROUND_ELEMENT", "North Arrow")[0]
mf = northArrow.mapFrame
northArrow.elementPositionX = mf.elementPositionX + 8.8
northArrow.elementPositionY = mf.elementPositionY + 0.7

# Align the title with the center of the map frame
title = lyt.listElements("TEXT_ELEMENT", "Name of Map Text")[0]
mf = lyt.listElements('MAPFRAME_ELEMENT',"Map Frame")[0]
title.elementPositionX = mf.elementPositionX + (mf.elementWidth / 3.7)
title.elementPositionY = mf.elementPositionY + (mf.elementHeight / 0.98)

# Reposition the Legend and fix legend title
legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
legend.title = "Legend"
legend.elementPositionX = mf.elementPositionX + 4.8
legend.elementPositionY = mf.elementPositionY + 6.15

# setting layout to bookmark
aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout
lyt = aprx.listLayouts("Layout_King")[0]
mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark
bkmks = mf.map.listBookmarks()
bkmks.sort(key=lambda x: x.name, reverse=True)
for bkmk in bkmks:
      mf.zoomToBookmark(bkmk)
      lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")

Is there something I'm doing wrong, or missing?

Also I'm using this custom script tool ModelBuilder

I tried using legend.updateItem() and legend.updateItem("*"), but it gives me the following error:

AttributeError: 'LegendElement' object has no attribute 'updateItem'


texts: What software would be useful to create General Transit Feed Specification (GTFS)?

texts: I am doing a project in which i am required to create GTFS. I look at transit google developer forum which provide a format to create feed but that completely have to done by hand (we have to write text files), i was looking for a software that can help in this matter.


texts: Query Geopackage using spatial index in QGIS DB manager

texts:
I would like to do spatial join in QGIS DB manager on geopackage database.
The query should look like this

select * from Soil_C_cadaster a,
(select * from Soil_C_cadaster
where forest_c is null) b
where st_touches(b.geom,a.geom)

It's a cross join on the same layer with 13000 not so complex features but I'm limiting it to some 200 features in table b. Query takes more than two minutes, so I'm suspecting it's not using spatial index by default as would PostGIS with spatial predicate.

Here is example of tested polygons:

In post like this, there's written about use but just with one layer and search bounding box, not with more layers. I can't imagine how to write it correctly.


texts:
I have installed the OSGeo4W web server tools on a Server 2008 machine, and I would like to serve MBTiles to Leaflet.  I have tried these:

Zverik's mbtiles-php

Infostream's mbtiles-php fork

The required PHP files and .htaccess were placed in a directory under \htdocs and Apache happily serves files from there.  However, I am unable to form any URL that retrieves tiles.  Firebug shows a '200 OK' on the tile request, but no content.

Given that the mbtiles database and the script reside in a directory called 'mbtiles', can anyone inspect the scripts and explain how to properly form the URL?  A proper L.tileLayer string for Leaflet would also be accepted.

Also, any other minimal configuration solution (that will run on the OSGeo4W install) that will serve my .mbtiles to Leaflet would be accepted.


texts: Shift in the data between CartoDB and stamen layer, Google Maps API

texts:
I have a geojson loaded into CartoDB database. And I have a CSV file too. Names field would be common in both the GeoJSON and CSV file. I had written a logic if the place name in CSV matches the place name in GeoJSON then just add the corresponding geometry from the GeoJSON to the stamen maps. Now my problem is the GeoJSON data matches perfectly after it is loaded into CartoDB. But when the geometry from GeoJSON is added to stamen maps (if the place name of CSV matches to that of GeoJSON ) there seems a slight shift in the data in stamen maps. I think this would be due to the extent had not been set properly.I would like to know how google maps api is setting extent at the backend for displaying data? What would be the exact problem? The snapshot of data inside CartoDB and stamen maps are attached below.

Data in CartoDB

Data in stamen maps

This is the code which am following,

 var data = {
                "type": "FeatureCollection",
                "features": feature
            };
            var layer = "toner-lite";
            var styles = [
                {
                    featureType: 'water',
                    elementType: 'geometry',
                    stylers: [
                        {hue: '#000000'},
                        {saturation: 100},
                        {lightness: -100}
                    ]
                }
            ];
            var mapProp = {
                center: new google.maps.LatLng(38.75597, -77.974228),
                zoom: 4,
                disableDefaultUI: true,
                zoomControl: true,
                zoomControlOptions: {
                    position: google.maps.ControlPosition.LEFT_TOP
                },
                mapTypeId: layer,
                mapTypeControlOptions: {
                    mapTypeIds: [google.maps.MapTypeId.ROADMAP, 'map_style']
                }
            };
    var styledMap = new google.maps.StyledMapType(styles,
                    {name: "Styled Map"});
            var map = new google.maps.Map(document.getElementById("googlemaps"), mapProp);
            var geocoder = new google.maps.Geocoder();
            var bounds = new google.maps.LatLngBounds();
            map.mapTypes.set('map_style', styledMap);
            map.setMapTypeId('map_style');
            //map.mapTypes.set(layer, new google.maps.StamenMapType(layer));

            map.data.addGeoJson(data);

And in the above code the object data has a value "feature" for the key "features". The structure of the value feature would be as below


texts: coordinate returned by xyCoordinate signal is different

texts:
I am working on a project based on the api qgis 2.8 and qt 4.8.5
I have two computer one with the qgis 2.8.1 and the other with the  qgis 2.8.3
In my application I have a connect to the signal xyCoordinate to get the mouse position on the map

connect( this,  SIGNAL(xyCoordinates(QgsPoint)),
 this,  SLOT(displayCoord(QgsPoint))    );

I have two different result
on the 2.8.1 the coordinate returned is corresponding the latitude and longitude of the point
but on the 2.8.3 the coordinate is corresponding to nothing it look like double number

the code of my application is the same on both computer, same system ubuntu 14.04lts, qt4.8.5 libs, the same map

is there a difference between the qgis 2.8.1 c++ api and the qgis 2.8.3 ?


texts: How to output a raster of special case's distribution with regard to overall distribution?

texts:
For example, the crime scenes are unevenly distributed in a city. Each case's scene can be represented as a point. Thus, we can generate a raster map, which shows the crime scenes' density in this city.

Also, we can draw such a type of raster for any particular crime type, like homicide.

Suppose we have generated a density raster for overall cases and a density raster for homicide only. Now I want to find the correlation between the distributions of homicide and overall cases. In other words, I want to find the place which have a significant higher homicide risk, with regard to overall crime distribution in this city.

Now I can only think about the raster calculator in ArcGIS, but I am afraid if it oversimplifies my issue. Is there any other tool in ArcToolbox could better solve my question? Any recommendation is welcomed.


texts:
I want to update the mshp-file of a shapefile with SAGA.
My issue is: I renamed some columns of a shapefile with SAGA (saga_cmd table_tools 23). For renaming I had to give SAGA the dbf-file, because it doesn't work with a shp-file. But if I now create buffers for this shapefile, the buffer has the old field names. If I delete the mshp-file, everything works. However, if possible I would like to keep the mshp-file.
Is there any way to update the mshp-file? When looking at the files in a text explorer, the dbf-file contains the new field names, but the mshp-file contains the old field names.


texts: Suggestions of route optimization algorithms to implement in QGIS

texts: I need to compute routes on a QGIS project with time and distance associated. Is there any route optimization algorithm to compute routes given the points that one wants to connect?


texts:
I was doing a query to extract pixel value from a raster and came up to some strange results. If I extracted the value using R and the Raster package, the value was different from the same point extracted using PostGIS.

Doing some further investigation, I found out that the differences only occurred when I used 'round' coordinates, that is -55.0, -22.0. These coordinates are in between two columns in my raster. So what I figured was that R was choosing the pixel from the right side of the edge and PostGIS, from the left.

Does each software round the coordinates in a different way?


texts: Copying linestring attribute column to start_point end_point attribute table using PyQGIS?

texts:
I'd like to extend the solution proposed to Can vector layer get start point and end point of line using PyQGIS? and retrieve linestring attributes (columns) to the start and end points from the original object.

In the example below the Polyline object have a column [ID_Rota] that I'd like to copy to the resulting point table.

I have - so far unsuccessfully - tried to implement @joseph's recommendation and add a feat.setAttributes(feature.attributes()) with:

line_layer = qgis.utils.iface.activeLayer()
feat = QgsFeature()

point_layer = QgsVectorLayer("Point?crs=epsg:4326", "point_layer", "memory")
pr = point_layer.dataProvider()

for feature in line_layer.getFeatures():
    geom = feature.geometry().asPolyline()
    start_point = QgsPoint(geom[0])
    end_point = QgsPoint(geom[-1])
    feat.setGeometry(QgsGeometry.fromPoint(start_point))
    pr.addFeatures([feat])
    #feat.setGeometry(QgsGeometry.fromPoint(end_point))
    #pr.addFeatures([feat])
    id_rota = feat.setAttributes(feature.attributes(ID_Rota))
    pr.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(point_layer)

and the error

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "~/temp/tmph65dmx.py", line 15, in <module>
    id_rota = feat.setAttributes(feature.attributes(ID_Rota))
NameError: name 'ID_Rota' is not defined


texts: Cutting a polygons into smaller squares in Python

texts:
In order to analize a GeoTIFF image we'd like to cut a huge polygon feature from a shapefile into smaller chunks (this is because for some reason our analysis is not performed well with a polygon with a huge extent).

What exactly we want to obtain is to for a given polygon P, window_width and window_height we want to get a list of polygon features where every feature is a intersection of a P with a set of disjoint, adjacent windows which has a shape of a rectangular with window_width and window_height.

This picture shows what we mean:

For this polygon we'd like to get a list of 12 new polygon features where first is a triangle, second - a pentagon, third - a square and so on.

What is important that we want to have this as a Python function. So we'd like to have the following function:

cut_polygon_into_windows(p, window_height, window_width)

to be add to our API.


texts:
I am trying to simplify a polygon layer of parliamentary constituencies in Kenya using SimplyPy in QGIS and appear to be getting an erroneous result.

The method is working for most of the polygons but there are 10 or so which are disappearing in the final result. Interestingly, the polygons which are disappearing all have similar object IDs (from 280-290 - there are 290 altogether) and so it could be due to SimpliPy having a limit of only 280 polygons at once... (I hope not!). Before and after below:

I have included a link to my shapefile here.

My question is: Why are my polygons disappearing in the final output?


texts: Merging features from equal attributes in QGIS 2.18.3?

texts:
I am working on a feature merge Operation, able to stick features together, if they have the same attributive values (e.g. TP_RKZ1 & TP_RKZ2 are empty so stick the multiple features to one single feature together if the following feature is not equal to the feature before [like in the Image]).

To get this behavior, I created the following code, which is actually already selecting the features in the mentioned way. The only missing part is the "merge process". Here I am not able to figure out the misbehavior of the used geom = geom.combine(feature.geometry) syntax.

Any suggestions?

"""APPROACH"""
"""Loop comparing attribute ID_RKZ1 & ID_RKZ2 values of the actual feature with the following feature and add the feature to the selectedFeatures using the
   myWorkLay.select () function and if it is not equal add last equal feature to the selectedFetures
   and perform the merge of all equal features"""
for i in range(0,len(myAttributes)-1):
    if myAttributes[i][1] == myAttributes[i+1][1] and myAttributes[i][2] == myAttributes[i+1][2]:
        print "My items {0} | {1} are equal".format(myAttributes[i],myAttributes[i+1])
        'Select feature because of feature-id'
        myWorkLay.select([myAttributes[i][0]])
        #print "%s"%myAttributes[i][0]
    else:
        print "!My items {0} | {1} are not equal!".format(myAttributes[i],myAttributes[i+1])
        'Select feature because of feature-id'
        myWorkLay.select([myAttributes[i][0]])
        print myWorkLay.selectedFeatureCount()
        #print myWorkLay.setSelectedFeatures()
        'Iterate through all selected features and merge them'
        for feature in myWorkLay.selectedFeatures():
             print feature.geometry()
             if geom == None:
                geom = feature.geometry()
             else:
                geom = geom.combine(feature.geometry())
        myWorkLay.removeSelection()
        geom = None

Here the result of my selection and

 10 ' => result of print myWorkLay.selectedFeatureCount()
 ' following the output of print feature.geometry()
<qgis._core.QgsGeometry object at 0x000000001C1742F0>
<qgis._core.QgsGeometry object at 0x000000001C174400>
<qgis._core.QgsGeometry object at 0x000000001C174378>
<qgis._core.QgsGeometry object at 0x000000001C1742F0>
<qgis._core.QgsGeometry object at 0x000000001C174400>
<qgis._core.QgsGeometry object at 0x000000001C174378>
<qgis._core.QgsGeometry object at 0x000000001C1742F0>
<qgis._core.QgsGeometry object at 0x000000001C174400>
<qgis._core.QgsGeometry object at 0x000000001C174378>
<qgis._core.QgsGeometry object at 0x000000001C1742F0>


texts: Importing GeoTIFF into PostGIS without raster2pgsql?

texts:
I know that one can import GeoTIFF into PostGIS via raster2pgsql. I was wondering if there is any easy way (as a sequence of SQL statements) to do the same without resorting to the command line tool?

Essentially, I am interested in a PostgreSQL function that is the inverse of ST_AsTIFF, which given a pathname to a GeoTIFF, creates a pgraster table or a row in such a table.

The reason I ask is that I need such a function in a SQL stored procedure to import and process GeoTIFFs (not known in advance).


texts:
I have .shp and .csv files for polygons that I have extracted from the georeferenced version of the pangea supercontinent. Link to the image which I have started from: http://4.bp.blogspot.com/-j6ZGqXskUm8/UDem8twbzaI/AAAAAAAACfY/fs-E3_PR-N8/s1600/pangea_politik.jpg

In mapbox: I am trying to use these files to make an underlying custom map, to be further used in Tableau. When I uploaded the .shp file as a tileset, it was saved as a vector layer. I imported it into the 'style', and tried to 'publish' it. The published map could not link with the data that I used in Tableau, and the heatmap that I published used the country shapefiles in the current day, and not the pangea ones that I made.

I want to finally produce a heatmap of users on the pangea map. Can somebody help me out?


texts: Custom JSONDecoder for ID Task Result, ESRI Flex API

texts:
When the ESRI API decodes a response from the identify task, there is a few second lock in the Flash application which is due to the decoding of the results.

I have used chunk parsing algorithms in action script to prevent the UI from locking, and I would like to implement this for the ID task response decoding.

I assume the class I need to extend is the com.esri.serlization.jason.JSONDecoder class. Does anyone know how I can get that code so I can re-write it, and use this extended class to parse the results from the ID task? Or, if I use a third party JSON decoder and re-write it, how can I use this class to parse the results from the ID task?


texts: How do I get the coordinates of my data frame always in long/lat?

texts:
I want to determine the coordinates of my current data frame. To do so I wrote this script:

def DataFrame ():

  mxd = arcpy.mapping.MapDocument("CURRENT")
  df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
  return mxd,df

def Center (df):

  extent = df.extent
  xmin = extent.XMin
  xmax = extent.XMax
  ymin = extent.YMin
  ymax = extent.YMax

  x = xmax - ((xmax - xmin)/2)
  y = ymax - ((ymax - ymin)/2)
  return x,y

As long as I am working with a geographic coordinate system, the results will be exactly the longitude and latitude. But as soon as I have a projected coordinate system it will return map units (i.e. meters). How can I avoid that problem? I always need to get lon/lat independend on the maps projection.


texts: Only intersect shapefiles that overlap arcpy

texts:
I have multiple shapefiles in a folder and I want to intersect all of them which overlap.  So far I am making two lists and overlapping all of the possible file combinations (I am using a cartesian list in this example, although in hindsight I think I could have skipped that step).

I am using this code:

import arcpy, os, itertools

arcpy.env.workspace=r'F:\Sheyenne\fire\split_fire'
shapefiles=[os.path.join(arcpy.env.workspace, i) for i in arcpy.ListFeatureClasses('*.shp')]
shapefiles2=[os.path.join(arcpy.env.workspace, i) for i in arcpy.ListFeatureClasses('*.shp')] #same shapefiles but used to find the file combinations
out=r'F:\Sheyenne\fire\intersects'

if not os.path.isdir(out):
    os.mkdir(out)

result = list(itertools.product(shapefiles,shapefiles2))

for files in result:
     x=(files[0])
     y=(files[1])

     for i in x:
         outname=x.split('\\')[-1].split('.')[0] + '_' + y.split('\\')[-1]
         arcpy.Intersect_analysis([x, y], os.path.join(out,outname))

print "Done"

This will intersect all the possible combinations, but since there are so many combinations it is extremely slow (the majority don't overlap), I am thinking there is a faster way to determine which ones instersect in the first place and then only instersect those ones.  My code above also doesn't deal with the fact that more than 2 shapefiles can overlap.


texts: What extensions are included in ArcGIS 10.2 Basic Version

texts:
I wanted to purchase ArcGIS 10.2. i used its trial version for last 2 month. now i am thinking to have Basic version. My work is based on Spatial analysis. so i dont need 3d and business analyst. for that i wanted know what extension will come with Arcgis 10.2 basic version.?. what are the functions available? and Spatial analyst extension come with ArcGis basic? or i have to purchase spatial analyst separately

if anyone have info about this plz tell me.


texts: Changing CRS of geometry in geopandas dataframe efficiently?

texts:
I have a shapefile which I am converting into a geopandas dataframe. The CRS of this shapefile is epsg:4326 and I want to convert this to epsg:3857.

import geopandas as gp

file_path = 'zip_poly.shp'
sf = gp.read_file(file_path)
sf = sf[['ZIP_CODE','PO_NAME','STATE','geometry']]
sf = sf.to_crs(epsg=3857)

You can find the shapefile here.

The shapefile has around 33K zipcodes and it's taking too long to convert. Is there a better way to convert quickly?


texts: Crop one shapefile using another shapefile in QGIS

texts:
In QGIS 2.0.1, is there a way to automatically slice up one shapefile and export it as multiple files using polygons defined in another shapefile?

I have a detailed shapefile that categorizes the land ownership of the entire state of Arizona.  (It's called "az_surf_mgmt" on this page.)  As a whole, it far exceeds the size and complexity restrictions of KML support in Google Maps.  So I want to take that file and chop it into individual counties, or perhaps game management units, as defined by another shapefile obtained elsewhere.

I found an almost identical question, but the accepted answer does not address how this might be done automatically.  If the answer requires that I learn Python and PyQGIS, so be it.


texts:
I run QField on both a tablet and my phone. In the last couple of days I suddenly can't open any of my projects on the tablet. They all work fine on my phone which is running the same version 1.6.4.  All of them share the majority of their files, which consist of a half dozen spatialite files and geopackages.
Last night I recopied several files to the tablet and opened the project I need today, but now that I am in the field it doesn't work. I've reinstalled the program as well.
To be specific this is what happens. Qfield starts normally and will open the sample projects that come bundled with it. If I select any of my projects it crashes during the project load. There are no error messages.
I cant find a log file in the ch.opengis.qfield directories. How can I find out what the problem is?


texts: How to cancel SaveEdits/StopEdits button press (ArcObjects)?

texts:
I have an interesting situation.Suppose a user creates 10 features in ArcMap in an edit session. Now the user tries to save this edits by clicking save edits or stop edits button. I would be performing some checks at this point(by listening to  OnStopEditing/OnSaveEdits events). If the result of this check is not satisfactory, I would want ArcMap not to commit these changes. Also the state of edits not should be the same as it was before the Save Edits/Stop Edits button was pressed.

So in the case of failed checks, the 10 digitized features would stay ther in map, without permanantly committing to the database.

I hope  I am not complicating things. Another way to describe the situation is, If the validations fail, ArcMap should just behave as if the user had not presed the StopEdits/SaveEdits button.

How can I do this using ArcObjects? I tried as shown below. But it simply saves edits and then stops editing. Anyone got ideas?

void editor_OnStopEditing(bool save)
{
  if (save == true)
  {
    if(HasMyValidationFailed())
    {
      //do not stop edit here..
      return; //did not work..
      //throw; //did not work..
    }
  }
}


texts: Toggle visibility of features in FeatureLayer

texts: How to toggle the visibility of features in a feature-layer using check-boxes. I have grid containing check-boxes in a columns representing each feature in Feature-layer and so if I check off the check-box I can turn off the visibility and if I check in turn on the visibility using ArcGIS JavaScript API 3.1.


texts: Why does buffer produce same error consistently?

texts:
I have created a point --> Figure 1

I have then selected the Buffer tool and created a buffer around my chosen point, with a degree interval of 0.012. I have selected the End style cap as square and run in background. This has produced a buffer around my point

Figure 2

Because I am using WGS84 CRS, I cannot use meters as the chosen map unit, instead I am using degrees and I have to guess a degree to match close to what I want. Does anyone have a solution to this specific problem?

My next problem is that after trial and error is that once created a buffer around the right dimensions they are constantly providing contradictory data.

Figure 3

From the right side it records a distance of 705m until the center, giving a diameter of around 1400m

Figure 4

From the left side it produces a distance to the center point of 1343m until the center and a diameter of around 2700m.

Figure 5

indicates that from corner to corner the distance is almost exactly 3000m.

Can someone please explain to me how buffers work. because I am using the Buffer tool the way people are telling me too, yet the results from it seem so inaccurate.


texts: Data about administrative subdivisions capital cities

texts:
I am looking for data about second or third rank administrative capitals in every African country (the administrative centre of a country territorial subdivision), but I can't find anything close to what I'm looking for.

So far, I found data on administrative boundaries (GADM, Natural Earth, UN SALB Project, OSM) and about cities (Natural Earth, OSM, Geonames), but nothing that specifically links an administrative subdivision to its administrative capital.

Although Geonames is supposed to provide that kind of information, I briefly checked and found that the data is of very poor quality.


texts: Spatial join to kernel density / kernel density with polygon/point output in Arc

texts:
I'm looking for a way to join polygons to a density function for some other feature of interest.

Current workflow

make density raster from feature set of interest, using point/line/kernel density
sum raster inside output geometry (for polygons) or just sample it (for points)

This seems inefficient as I have to calculate a density raster for many points I won't even use.  Is there a tool that does the same thing but without the raster?  This would be like a kernel-density-aware spatial join, or a kernel density function that outputs to a preexisting geometry rather than to a raster.


texts: Finding overlay/intersect shape percentages using QGIS?

texts:
I am working on a project in which it's necessary for me to overlay 2 layers (a map of precincts and a map of census tracts).

I need to be able to determine the %area of each census tract covered by each precinct. For example, if Census Tract 1 is 25% Precinct 1, 40% Precinct 2, and 35% Precinct 3.

How do I go about doing this?

I am using QGIS 3.8 and have the necessary shapefiles.


texts: Identifying when value in field changes from 0 to 1 in ArcGIS Desktop?

texts:
I am trying to identify the start and stop of a particular incident- these are indicated by the unique values 0,1,2,3,4,5,6.

0 represents no incident- therefore in order to identify when an incident starts and stops and only retrieve these data records (i.e. incident 1) I would like to know how to query or determine when the value in the incident field changes from a 0 to a 1 (i.e. the start) and query when the value changes from a 1 to a 0 i.e. the finish. I only want to compare these changes between adjacent records i.e. those either side of a record.


texts: Fine-tune label placement for road shields in QGIS

texts: I want to display road shields and road numbers within them. I already checked the option to connect multiple lines in order to avoid double placement, set the minimal length for lines to be labeled to 5mm and I even limited the absolute number of sign to 25 but it still looks cluttered on some roads whereas some have no sign at all. Even if restrict the above mentioned settings even more, I still get weird looking double labels like in picture 2.


texts: Spatial Join works when invoked through ArcMap GUI, but "Out of Memory Error" when using ArcCatalog "Spatial Join" tool?

texts:
(This question was inspired through the discussion here Spatial Join works when invoked through ArcMap GUI, but "Out of Memory Error" through ArcPy)

I want to join an input feature class (~11,000 polygons) with another feature class (~31,000,000 polygons). When I do this trough ArcMaps 'Join - based on spatial location'  tool, it takes quite long but eventually succeeds (adding 'Sum' for all joined fields). I monitored memory consumption throughout and it seems 'healthy' to me (not even going close to my 64 GB RAM).

When I try to do the same with the 'Spatial Join' tool:

I run into an "Out of Memory" error after a few minutes. What's strange about it is that, again, as far I can tell from the task manager, it's nowhere close from actually running out of my memory.

So, I believe I'm doing exactly the same thing through ArcCatalog as through the ArcMap GUI, yet getting no result. What is the analog of the "Join Data" in terms of ArcCatalog functions?


texts:
I need to make a heatmap with a bunch of linestrings using AngularJS and Open Layers 3. I didn't find any examples, just a question answered here but that don't worked for me. In the question he do a action listener to when the features are added to the layer but for me, I already have all the features so I don't need this action. I've made some alterations to fit in my project structure, he don't show any errors but he don't show the heatmap either!
I have a service that creates all the layers on my map, I need this heatmap to be a layer, so that way I can set it visible or not.
LayerService:

/**
 * Routes is the list of routes that contains the linestring
 * inside the "geom" property.
 */
LayerService.routesHeatMap = function (routes, blur, radius) {
  const centroids = [];
  routes.forEach(route => {
    const geometry = route.geom;
    if (!!geometry) {
      const point_feature = new ol.Feature({ weight: route.expansion });
      point_feature.setGeometry(new ol.geom.LineString(geometry.coordinates));
      point_feature.setProperties({ projection: "EPSG:4326" });
      centroids.push(point_feature);
    }
  });
  let vector = new ol.layer.Heatmap({
    title: 'Routes. Heatmap',
    source: new ol.source.Vector({
      features: new ol.Collection(centroids),
      overlaps: true,
      opacity: 1,
      visible: false,
      format: new ol.format.GeoJSON({
        dataProjection: 'EPSG:4326'
      })
    }),
    zIndex: 1,
    gradient: ['#007700', '#6ce200', '#efef00', '#ffff00', '#ff6600', '#ff0000'],
    blur: parseInt(blur.value, 10),
    radius: parseInt(radius.value, 10)
  });

  return vector;
};

My controller that have my map:

var blur = document.getElementById('blur');
var radius = document.getElementById('radius');
heat = LayerService.routesHeatMap(response, blur, radius);
vm.map.addLayer(heat);


texts: Generalize point data to a reference vector grid

texts:
I have a point feature, each point has a floating value assigned. I have a square grid (polygon) feature overlapping the points. I want to assign a value to each grid cell equal to the sum of the points' value inside the cell.
Something similar to Zonal Statistics tool (ArcGIS) that can use a point layer instead of a raster for input data.


texts: Reprojecting Natural Earth raster into natearth projection with gdalwarp causing black squares?

texts:
I'm trying to reproject the Natural Earth 2 raster layer (http://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-2/) into the natural earth projection (http://en.wikipedia.org/wiki/Natural_Earth_projection) using gdalwarp.

The whole thing looks pretty good except I get 4 black squares by the edges of the map as seen below.

gdalwarp options used [gdalwarp included in QGIS: GDAL 1.11.0, released 2014/04/16]:

gdalwarp -overwrite -s_srs EPSG:4326 -t_srs "+proj=natearth +wktext" -r lanczos -of GTiff NE2_HR_LC_SR_W_DR.tif reprojected.tif


texts: Error while running Variable distance buffer in QGIS 3.12

texts:
I have problem in QGIS 3.12 when I try to use the tool Variable distance buffer I get an Error.
This is the log message I get from the error:

  { 'BUFFER' : 'TEMPORARY_OUTPUT', 'DARC' : 5, 'DISSOLVE ' : False,
  'DIST_FIELD' : 'raduis', 'DIST_SCALE' : 1, 'NZONES' : 1, 'POLY_INNER '
  : False, 'SHAPES' : 'C:\Users\Public\shape5.shp|layername=shape5' }

  shapes_tools "Shapes Buffer" -SHAPES "C:\Users\Public\shape5.shp"
  -DIST_FIELD "raduis" -DIST_SCALE 1.0 -NZONES 1 -DARC 5.0 -DISSOLVE false -POLY_INNER false -BUFFER
  "C:/Users/לירון/AppData/Local/Temp/processing_sFcGXG/44994fdb08b0468883ff9936f0c2a56e/BUFFER.shp"

  Execution completed in 10.89 seconds Results: {'BUFFER':
  'C:/Users/לירון/AppData/Local/Temp/processing_sFcGXG/44994fdb08b0468883ff9936f0c2a56e/BUFFER.shp'}

How can I solve it?


texts:
Looking at OpenStreetMap's Transport layer, I can see public transit routes, but not any ticket zones (for example Berlin's ABC zones or London's 1-9 Underground zones).

Does this information exist on OSM, and if not are there any proposals for specific tags to add it?


texts: Assign attributes from nearby line to another line in FME

texts:
I have two polyline datasets. One is a topoligically/spatially correct road centerline (red in the image below "new"), the other is a spatially inaccurate but properly attributed road centerline (yellow - "old"). I want to get the attributes into the spatially accurate one.

I am just learning FME and would like to do it using that.

It's basically a spatial join but with a tolerance of up to about 5m . I see that there are options for point datasets (Spatial join using FME) but I can't seem to find anything for polylines.

I found a solution based on arcmap but I prefer to use FME in this case. See Merging and joining two line feature classes while maintaining attributes from both using ArcGIS for Desktop?

There's the following but it also dealing with a point dataset (FME - Is there a transformer to snap point on point of nearest line?)

Do I need to buffer the new dataset and then join the attributes of anything within the buffer and then convert this back into a line? or do I need to convert line to point and then run the above-suggested processes? or is there some simpler fme model process I can apply?

====
Based on @buddiebubba response I had the issue that he and @Mark Ireland said would happen.

I am not sure what attributes to use in the following

=== with the Azimuth Calc I don't get any matches...

See https://drive.google.com/file/d/0B1XGDTNNfbMJZXZOZWRKa2lwMG8/view?usp=sharing for the fmw

=== I also tried the following which ends up with pretty much the same result.

See https://drive.google.com/file/d/0B1XGDTNNfbMJXzBKaDNWZldxNGM/view?usp=sharing for the fmw

=== the fmw works but there's still a lot of mismatch...

The sample data is in https://drive.google.com/open?id=0B1XGDTNNfbMJVzI1TUdEOFg2Yzg


texts: How to color only the borders where the polygons are adjacent?

texts:
Is there a way to only color the borders where the polygons are adjacent?
For example, I want to color the borders in the green circle, but not the ones in the red cirlce.


texts: Matching (joining) two slightly incompatible polygon shapefiles

texts:
I have two shapefiles that represent two slightly different zonings of the same area (municipality zoning and transport zones), about 3000 polygons each. Several smaller municipalities may form a single transport zone, while a larger municipality may be split into several transport zones.

To these shapefiles, two datasets are linked: Traffic flow is available for the municipality zoning, and travel time is available between transport zones. I need to approximate the travel time between municipalities. My idea is to estimate the travel time between two municipalities as the weighted average of the travel times between all related transport zones. The weighting would happen according to the share of area that a transport zone contributes to the corresponding municipality.

I have come up with an R code that seems to do the right thing. However, I've also read about the MAUP (Modifiable Areal Unit Problem), and it seems to apply here. My questions are:

Is the construction of a spatial mapping between incompatible zonings a common practical problem?
If yes, how is this problem solved in practice -- in R, QGIS, ArcGIS (How to merge two 'incompatible' polygon layers in ArcGIS?) or using other tools?
Is there a more sensible weighting, perhaps one that incorporates centrality (i.e., those transport zones that are located more centrally would have more influence on the result)?

References to academic publications are more than welcome.


texts: Adding multiple cartodb layers in viz.json to Leaflet layer control?

texts:
I have create a visualisation in cartodb with 4 layers. The resulting viz.json is added to a leaflet map as follows:

cartodb.createLayer(cartomap, [url], {https:true})
                .addTo(cartomap)
                .on('done', function(layer) {
                // add individual layers
}

I want to accomplish that all 4 layers in the viz.json can be toggled separatly in the layer control (switching both the markers/features and the legend)

The layer in on('done') is a valid leaflet layer and can be added to the Layer Control. However, it contains all sublayers which can then not be toggled individually.

The layer also contains an array with all the "sub" layers in there, but these are not valid Leaflet layers (onAdd function is missing), so these can not be added to the layerControl individually. The layer.getSubLayer(i) also returns a not valid Leaflet layer.


texts: Using Append tool in ArcPy with Edit Session?

texts:
I need to run the append tool in an edit session because the feature class I am appending to is in a gdb with a relationship class, and I need to update the field mapping for a few of the fields. I created my Append tool using ModelBuilder to get the "field_mapping" correct.

I created this code and it seems to work okay, however, when I incorporate it into my longer field (with this being the last step) it no longer updates the fields... However, no error seems to occur either.

Import arcpy module here
import arcpy
arcpy.env.overwriteOutput = True

# Script Arguments
Fitting = arcpy.GetParameterAsText(0)
Final_Layer = arcpy.GetParameterAsText(1)
Out_Workspace = arcpy.GetParameterAsText(2)

edit = arcpy.da.Editor(Out_Workspace)
edit.startEditing(False, True)

# Process: Append
arcpy.Append_management(Final_Layer, Fitting, "NO_TEST", "ANCILLARYROLE \"AncillaryRole\" true true false 2 Short 0 0 ,First,#;ENABLED \"Enabled\" true true true 2 Short 0 0 ,First,#;EDITEDBY \"Edited By\" true true false 64 Text 0 0 ,First,#,Final_Layer,W_Valve_EDITEDBY,-1,-1;DATEEDITED \"Date Edited\" true true false 8 Date 0 0 ,First,#,Final_Layer,W_Valve_DATEEDITED,-1,-1;PLACEDBY \"Placed By\" true true false 64 Text 0 0 ,First,#,Final_Layer,W_Valve_PLACEDBY,-1,-1;DATEPLACED \"Date Placed\" true true false 8 Date 0 0 ,First,#,Final_Layer,W_Valve_DATEPLACED,-1,-1;LABELTEXT \"Label Text\" true true false 64 Text 0 0 ,First,#;WATERTRACEWEIGHT \"Water Trace Weight\" true true false 4 Long 0 0 ,First,#;SYMBOLROTATION \"Symbol Rotation\" true true false 8 Double 0 0 ,First,#;IPID \"IPID\" true true false 4 Long 0 0 ,First,#,Final_Layer,Count_One_IPID,-1,-1;SUBTYPECD \"Subtype Code\" true true false 4 Long 0 0 ,First,#,Final_Layer,Count_One_SUBTYPECD,-1,-1;NOMINALSIZE \"Nominal Size\" true true false 8 Double 0 0 ,First,#;MATERIAL \"Material\" true true false 10 Text 0 0 ,First,#;BRANCHSIZE \"Branch Size\" true true false 8 Double 0 0 ,First,#;DEGREES \"Degrees\" true true false 8 Double 0 0 ,First,#;FRAMME_BASENAME \"FRAMME_BASENAME (FRAMME Tag)\" true true false 20 Text 0 0 ,First,#;FRAMME_IPID \"FRAMME_IPID (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_PRMRY \"FRAMME_RB_PRMRY (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_SCNDRY \"FRAMME_RB_SCNDRY (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_FSC \"FRAMME_RB_FSC (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_STATE \"FRAMME_RB_STATE (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_COMPONENT \"FRAMME_RB_COMPONENT (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;FRAMME_RB_OCCURRENCE \"FRAMME_RB_OCCURRENCE (FRAMME Tag)\" true true false 8 Double 0 0 ,First,#;GPSX \"GPS X Coordinate\" true true false 8 Double 0 0 ,First,#;GPSY \"GPS Y Coordinate\" true true false 8 Double 0 0 ,First,#;GPSYESNO \"GPS Yes/ No\" true true false 255 Text 0 0 ,First,#", "Plug")

edit.stopEditing(True)

The only difference between this snippet and my larger script, is that my larger script has different parameters, the "Final_Layer" is actually created in the script rather than being an input... But I don't think that would make much of a difference.

Is there a reason the "editing" session is not kicking in during m larger model?

Or is it related to what I have here?


texts: Concise learning resourses for OGC Filters or Filter Encoding Standard (fes)

texts:
When I am doing HTTP GET requests to GeoServer it is fairly easy, but when I need to do http POST request I need to pass XML formatted parameters in the body and for filtering I need OGC filters and filter encoding standard of OGC.

So I want to learn it, so from where should I start and what are the resources for that kindly exist?


texts: How do I write a script to list those layers with broken data sources?

texts:
In my SHP to GDB toolbox, I create a list of layers in MXD, then I check for those not broken I move them. Despite using ListBrokenDataSources I have errors like this:

  Runtime error : [Error 3] Impossible
  trovare il percorso specificato:
  u'J:\20120208_Backup_disco_E\Archivio\ATLArete\Shape Atlarete
  1-1-2008\miste.shp'

Translate:

  Runtime error : [Error 3] Unable to
  find the specified path:
  u'J:\20120208_Backup_disco_E\Archivio\ATLArete\Shape Atlarete
  1-1-2008\miste.shp'

Here is my code:

# Create a list of SHP in mxd
    list = []
    df_list = arcpy.mapping.ListDataFrames (mxd)

    for df in df_list:
        for fc in arcpy.mapping.ListLayers (mxd, "", df):
            if (fc.supports("DATASOURCE")) and (fc.dataSource.endswith(".shp")):
                        list.append(fc)
    print("List created")

    # Check for borken datasource
    for fc in list:
        if arcpy.mapping.ListBrokenDataSources(fc):
            list.remove(fc)
    print ("Broken data removed from list")

    # Move fc to gdb
    add_list = []
    for fc in list:
        fc_date = time.strftime("%Y%m%d", time.gmtime(os.path.getctime(fc.dataSource)))
        try:
            dsc = arcpy.Describe(fc)
            out_temp = '_' + fc_date + '_' + fc.name + '_' + fc_prj
            out_name = arcpy.ValidateTableName(out_temp)
            print ('Created ' + out_name + ' to ' + gdb_full_path)
        except:
            out_temp = '_' + fc_date + '_' + fc.name + '_' + 'UNKNOWN'
            out_name = arcpy.ValidateTableName(out_temp)
            print ('Created ' + out_name + ' to ' + gdb_full_path)
        arcpy.FeatureClassToFeatureClass_conversion(fc, fd_path, out_name)
        replace_name = fd_path + '\\' + out_name
        fc.replaceDataSource (gdb_full_path, "FILEGDB_WORKSPACE", out_name)
        print ('Replaced data source ' + out_name + ' in ' + gdb_full_path)

    # Create a new MXD
    new_mxd = os.path.splitext(gdb_full_path)[0] + '\\' + 'MXD' + '\\' +  time.strftime("%Y%m%d", mxd_date )+ '_' + (mxd_name)
    mxd.saveACopy(new_mxd)
    print ('Created new MXD ' + new_mxd + ' in ' + ws)


texts: How accurate are Ordnance Survey MasterMap topography spot heights?

texts:
The spot heights in the Ordnance Survey MasterMap sample data over Exeter with feature code 15404 are presented as an integer value in metres. But how accurate are these spot heights? I had assumed that spot heights are GPS measurements that are accurate to a few centimetres.

(First posted on the Ordnance Survey forums)


texts: Repairing a broken data source for a single layer in multiple map documents

texts:
I am a Python newbie and I'm trying to repair the data source for an aerial raster in close to 50 individual map documents.

My supervisor deleted the 2014 recently and uploaded the 2015 in its place. The new raster is in the same workspace and has virtually the same path, just a different name.

This is the template I was working with from ArcGIS resources. How would I need to modify it for the task at hand?

import arcpy
mxd = arcpy.mapping.MapDocument("C:\Project\Project.mxd") #I'd like to create a list

for lyr in arcpy.mapping.ListLayers(mxd):
      if lyr.dataSource == (r"C:\Project\Data\Aerial\SID_2014"):
        lyr.replaceDataSource(r"C:\Data\Aerial\SID_2015", "SHAPEFILE_WORKSPACE", "SID_2015")


texts: Determining which GeoServer version is running?

texts: Is there a way to find out which version of GeoServer is running from any standard document such as GetCapabilities?


texts: Error in calculating covariance matrix and fitting variogram model

texts:
When constructing variogram model, usually the spatial points are given, and we build a variogram model and conduct parameter estimation. But here in my question the points are not given while parameters for the variogram model is given, the spatial process is defined to be conforming to normal distribution with mean=0, but its sigma matrix is undefined and I have to calculate it using output values of variogram model. I calculated the gamma matrix using code:

h.Matrix<-as.matrix(dist(h.Matrix))
gamma.Matrix<-1-1/{exp(3*h.Matrix)}

But I need the value of C(0)(i.e.the value of sill) to caculate C(h) from C(h)=C(0)-γ(h). In this case, I suppose the value of sill is not defined?
And I have an error when fitting the variogram model as follows:

sp<-data.frame(x,y,z)
coordinates(sp)<-~x+y
v<-variogram(z~1,sp)
plot(v)
fit.wls <- fit.variogram(v, vgm(c("Lin", "Exp", "Sph", "Gau")))

Error in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges = fit.ranges,  :
fit.method 7 will not work with zero distance semivariances; use another fit.method value

I tried to change the fit.method but another error occured:

In fit.variogram(object, x, fit.sills = fit.sills, fit.ranges = fit.ranges,  :
No convergence after 200 iterations: try different initial values?

I cannot figure out what's wrong in the codes and would like to seek for your help...


texts: What is the official term for Photo Index?

texts:
I have a point shapefile that has all the aerial photos that were flown on our forest. However, they are very old and are black and white. They were scanned by our contractor. Their years are dated back to the 1930's to around 1950's.

I am not really familiar with the photo index and I have been searching around on the internet as I want to learn more about them. I want to know what they are for and how it works.

I have started working on them with my ERDAS Imagine 2016's AutoSync.

My question to you is when I google it, and  I have entered the word as  flight line or Photo index and it returned few results. Is there a official word/term that should I use to search for more information on google ?


texts: How to configure WFS layer with OpenLayers in PHP WampServer?

texts:
I have a problem to display WFS layer outside the GeoServer. My code working fine in GeoServer, but when I tried my code in WampServer(PHP), only WMS layer is being to displayed.

I've gone through many tutorial and stack overflow queries related to adding WFS layer in OpenLayers but nothing working in my case... I have also tried proxy.cgi scripts to allow proxy call to GeoServer, but it not works. Is CGI proxy script can be used only with Python or it can be used with PHP as well?


texts: Split multiple attributes from a big polygon to smaller ones

texts:
I am looking for a sort of spatial (un)join function.

Let's say I have 2 tables, both with a large number of polygons. The first one has a pop/polygon column in his attribute table. The 2nd one has smaller polygons that fit exactly into the bigger ones has a sum (these polygons represent census limitations).

Is there a function under QGIS or PostGreSQL that would separate the population/polygon from the 1st table into the smaller polygons of the 2nd table.

Example of the situation and end result: ID1 of table_1 has a population of 4. ID1, ID2, ID3 and ID4 of table_2 fit into ID1 of table_1. The function would create a population column into table_2 where ID1, 2, 3 and 4 each have a population of 1.

In other words I am looking for a function to do the exact opposite of the "Join attributes by location" in QGIS.


texts: arcpy.mapping.UpdateLayer returns runtime error

texts:
I am trying to update the style for a layer present in my saved map document with a source .lyr file present on my disk in the home folder. It returns an assertion error. The code I am trying to run is as follows:

import os, csv, arcpy
arcpy.env.workspace = "\ArcGIS\Python_test"
PointFile = "\ArcGIS\Python_test\Well_data.csv"
arcpy.MakeXYEventLayer_management("Well_data_View, "Easting", "Northing", "Well_plot_EventLayer", "WGS 1984 UTM Zone 43N")
arcpy.Delete_management("Well_data_View")
arcpy.FeatureClassToFeatureClass_conversion("Well_plot_EventLayer", "\ArcGIS\Python_test", 'Well_plot')
arcpy.Delete_management("Well_plot_EventLayer")
mxd.saveACopy("\ArcGIS\Python_test\Python_workspace.mxd")
mxd = arcpy.mapping.MapDocument("\ArcGIS\Python_test\Python_workspace.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
ul = arcpy.mapping.ListLayers(mxd, "", df)[0]
sr = arcpy.mapping.Layer("\ArcGIS\Python_test\Well_SR.lyr")
arcpy.mapping.UpdateLayer(df, ul, sr, False)

'PointFile' is a CSV table on my disk which contains XY information being used to plot 'Well_plot_EventLayer' which in turn is converted to a shapefile 'Well_plot' which I do not export to my disk. After saving my map document, I store the map document path to variable 'mxd', the map layer to be update to variable 'ul' and the source layer (style template) 'Well_SR.lyr' to variable 'sr'. The error generated by the last line of code is:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\mapping.py", line 1889, in UpdateLayer
    assert isinstance(source_layer, Layer)
AssertionError

Can someone please tell me where I've gone wrong?


texts:
I am using Qgis 2.1 ver. 5f5cd4c (64 bit version on win7 64) OSGeo installed. I was trying to reclassify a DEM in EPSG:32614. I used the SAGA's Reclassify algorithm in the processing toolbox. I defined that I wanted to use a simple table, defined the table and that I wanted the reclassified grid saved to a tif file. I also wanted the resulted reclassified DEM to be opened after the algorithm finished the calculations. This is what I am repeatedly getting:

I tried to open a new project and load the resulting raster manually, but Qgis does not recognize it as a valid raster data source:

I was able to open the file and convert it to GeoTIFF in SAGA. I wonder if I am missing something in the reclassifying procedure, or something related with the format of the resulting raster. Thanks in advance for any advice.

PS. I was able to use r.reclass to do what I wanted. I also tested the saga reclasify algorithm in other two installations and I have the same problem.

Gerardo


texts: Converting grid co-ordinates to GPS co-ordinates

texts:
I need to convert grid co-ordinates from an old document to GPS co-ordinates that I can plot into Google Earth and ArcGIS. Here are the instructions:

'It was unfortunately not possible to produce maps using UTM co-ordinates. This problem was overcome by imposing a grid onto the map of Zimbabwe and digitising the site maps. The origin (00.00,00) of the grid was at 25 degrees, 05 minutes, 24 second east and 22 30 00 S. A point (80.00,00) along the X-axis was at 32 50 40 E and 22 25 40 S. One unit of distance on the grid represented 10 kilometres on real ground. The co-ordinates given in the following table represent the positions of the sites within this imposed grid.'

An example of a grid co-ordinate is:
Name 2030 BD 1
Grid 60.21 24.78

Can anyone instruct me as to how I would go about doing this?


texts: What is a good method for serving many static map layers through ArcGIS Server 10?

texts:
I'm writing a Silverlight app that consumes some map services from ArcGIS Server 10.0 SP2 and I have about 15 layers of historical imagery I would like to give the user access to.  This imagery doesn't change and I'd prefer to use a cached map service so the user gets a smooth panning/loading experience.

The catch is the user needs to be able to have any two layers overlaid on each other so they can blend between the two.

I can see two options to do this:

Create 15 cached map services, one for each layer
Create one multi-layer cached map service

I've read several threads where it has been stated that multi-layer caches are not for web applications (and it appears they use the DynamicMapService from the API so they may not even satisfy all my requirements).

I'm sure this isn't an unusual situation to be in and I'm wondering what the typical solution might be.  Right now the layers are dynamic layers and it doesn't provide as smooth a user experience as I'd like.


texts: SQL-query in arcpy - SELECT WHERE feature between two values

texts:
I need to select features between two values in python with arcpy. However, whenever I try to run the code below, it returns with error message "An invalid SQL statement was used".

whereclause = """"X" BETWEEN {0} AND {1}""".format(
                                                  min(X)[0]-1000,
                                                  max(X)[0]+1000)
ap.FeatureClassToFeatureClass_conversion(SamledeMaalinger,ap.env.workspace,"SamledeMaalingerSort", whereclause)

The whereclause returns "X" BETWEEN 555157.0723 AND 557157.0723 so it's not the string that's the issue.

The SQL-query works if I use < and > instead, but the BETWEEN for some reason doesn't.

So what's the issue? Do the SQL-queries in arcpy have limitations?


texts:
I now need to reclassify, I'm using rastercalc... but I can't get the syntax right!!

I need that values >0 become 1

I found this example...
eq( [relief]@1, [mask]@4, 150 )

so I tried:
( [raster A insularis_1]@1, > 0, 1) and I absolutely failed


texts: Extract data within a range using ALOS DEM in GEE

texts:
I am trying to extract the LULC falling below 500m in Google Earth Engine. But not sure how we do that code here. Anyone who knows please help on this.
Using ALOS DEM here.


texts: Extracting qualified image patches in Google Earth Engine?

texts:
I want to extract qualified (e.g., cloud-free) image patches (e.g. size of 128x128) in a large image (e.g. 10000x10000) in Google Earth Engine. These image patches can overlap with other qualified patches. I will then feed these image patches to TensorFlow. Searching Earth Engine APIs, so far I only found ee.Image.clip, ee.Image.clipToBoundsAndScale, ee.ImageCollection.getRegion that may related to my goal. But they all require some geometry or feature or feature collections in latitude/longitude coordinates as the boundary parameter.

Is there a way to extract image patches in Earth Engine like slicing sub-array of a numpy array using indices (e.g., image[0:128, 0:128])?

Nicholas Clinton's overall flow works. I made some minor changes to the mask and the sample steps.

For the mask step, I used the ee.Image.fastDistanceTransform to calculate the Chebyshev distance to pixels with non-zero value (swap 0 and 1 in a normal mask, i.e., SwappedMask) and create the buffered mask by setting a threshold:

var Chebyshev_dist_to_Mask = SwappedMask
  .fastDistanceTransform({
    neighborhood: 46,
    units: 'pixels',
    metric: 'chebyshev'
  });
var BufferedMask = Chebyshev_dist_to_Mask
  .gt(ee.Image(33))
  .select(0)
  .toByte()
  .rename('buffered_mask');

For the sample step, I used ee.Image.stratifiedSample instead of ee.Image.sample to sample on good pixel class only.


texts:
Why are my GeoJSON files being saved without line breaks? Why are they missing the proper formatting? All features in the file are being saved in the same line.

I'm using QGIS 3.6.3.

Here is what I believe is "proper formatting"


texts: Colouring three individual attributes in point layer using QGIS?

texts: I am new to QGIS and am learning the expressions function. I would like to know how to independently colour three separate attributes on a point layer as below:


texts:
Is it possible to organize some shapes into folders to turn on/off all the shape togheter?
Is it useful for the shape that you can punt at the top/bottom of a project.


texts: writeOGR error: length(dsn) == 1L is not TRUE

texts:
I'm running a for loop R that results in a shapefile being written with each iteration. The for loop runs successfully once, with a shapefile being written out, but then on the second iteration I get this error:

Error in writeOGR(temp, outfile, layer = complete_paste, driver = "ESRI Shapefile",  :
  length(dsn) == 1L is not TRUE

For some reason, I cannot find any other question encountering this error with writeOGR(). I don't know what it means. Can anyone explain?

If it helps, here's the entire code I'm running. Basically I'm taking a large shapefile, and exporting smaller shapefiles from it by values in shp1$Vector.

library(rgdal)

shp1<- readOGR("/Users/JohnDoe/Desktop/Zone_Fixup/Z1/z1_scrub/z1_merge.shp")

output_path1<- "/Users/JohnDoe/Desktop/Zone_Fixup/Z1/z1_split/"

#Begin for loop to split shapefile by shp1$Vector values
for(i in 1:length(shp1$Vector)){

  #subset by join_id (Vector)
  temp<- shp1[shp1$Vector==shp1$Vector[i],]
  id_paste<- paste(temp$Vector)

  #Create strings for shapefile write-out
  outfile<- paste(output_path1, "Zone1_", id_paste, ".shp", sep = "")
  complete_paste<- paste("Zone1_", id_paste, sep = "")

  #Write out shapefile
  writeOGR(temp, outfile, layer = complete_paste, driver = "ESRI Shapefile", overwrite_layer = TRUE)

}


texts: CO column number density being large negative values (Google Earth Engine)

texts: I am trying to use Sentinel-5P NRTI CO: Near Real-Time Carbon Monoxide data for my project and this data is L3 processed data. The band description says that the estimated minimum value could be -279 mol/m^2 (the maximum value being 4.64 mol/m^2). I tried to find how the minimum value could be such a large negative value from the available documentation. But I couldn't find any source that explains the reason. Could someone explain?


texts: Lines appearing on vector layer after panning in QGIS

texts: One of my vector layers is creating random-seeming lines whenever I pan across the map. It's only one layer and the lines disappear if I turn it off. It doesn't appear in the print composer so it isn't ultimately that big of a deal, but it makes me wonder if something else is going on. Pictures are below, taken after panning just a bit between each one. Final pic shows the fill color blurs along with the border.