texts: AutoCAD Map 3D GeoTIFF in the wrong place

texts:
I recently was looking at a way of using GeoTIFF in AutoCAD 2015. I was looking at a way of specifying coordinate system using one of UTM projections. I didn't have much luck with that and I decided to try AutoCAD Map 3D 2017, there is an entire library of different coordinate systems including WGS84 UTM Zone 30N, the one I was interested in.

Now, let me explain my problem, I've got .tfw file that I generated using listgeo tool, here is the content of that file:

        4.9995042276
        0.0000000000
        0.0000000000
       -5.0003231419
    -9782.7177138953
  6572004.2399761314

I tried to use my GeoTIFF with ArcMap, QGIS and one more survey positioning package, all of them after specifying the right coordinate system place the GeoTIFF in the right geographical location apart from AutoCAD Map 3D.

When I specified the coordinate system and units in AutoCAD Map 3D, the entire GeoTIFF looks like it has been shifted a significant distance to the West.

I'm just wondering if AutoCAD Map 3D is drawing GeoTIFFs in a slightly different way or am I doing something completely wrong here?

Did anyone have any luck with using GeoTIFFs in Autodesk packages (Map 3d)? I'm worried that the .tfw file might have wrong values as in QGIS and that survey positioning package I didn't have to use that file and the software positioned GeoTIFF for me, when in AutoCAD Map 3D I actually have to place that file in the same directory as my GeoTIFF in order to be used but the x,y translation doesn't seem to be correct.
Here is the output of gdalinfo:

Size is 3805, 3248
Coordinate System is:
PROJCS["Transverse Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-2.716666666666667],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-9785.217466009140500,6572006.740137700900000)
Pixel Size = (4.999504227592880,-5.000323141986690)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=127
  TIFFTAG_YRESOLUTION=127
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   -9785.217, 6572006.740) (  2d53'17.62"W, 59d15'47.69"N)
Lower Left  (   -9785.217, 6555765.691) (  2d53'15.00"W, 59d 7' 2.84"N)
Upper Right (    9237.896, 6572006.740) (  2d33'16.92"W, 59d15'47.73"N)
Lower Right (    9237.896, 6555765.691) (  2d33'19.40"W, 59d 7' 2.88"N)
Center      (    -273.661, 6563886.215) (  2d43'17.24"W, 59d11'25.67"N)

@Jens
When I used my file in QGIS I didn't use any extra world files, just the geotif itself, however when I tried to load my file in AutoCAD Map 3D I had to use tfw file as without it I was asked to specify the offsets for my geotif image.


texts: Why does QGIS not connect to my new local PostGIS database?

texts:
I'm trying to connect to a new local PostGIS instance (PostgreSQL.app on mac) with QGIS 2.0, and it gives general connection errors.

Connection with same parameters with other clients (psql, PGAdmin) is ok, with both socket and port, so it can't be server configuration?


texts: Does Metadata for ArcGIS Python script tools in custom toolbox work?

texts:
I've created several python scripts and created tools from them.
I have a custom toolbox that contains all of my tools.
I add this toolbox to the ArcToolbox window in ArcCatalog.
I want to create metadata and set parameter explanations for each
tool.

I've noticed that, when I navigate to my toolbox using the ArcToolbox window, I cannot edit metadata by highlighting the tool in the toolbox window and choosing the 'Description' tab in the catalog window. I can, however, right click in the tool and choose 'Item Description', which opens another window and I can edit metadata in there.

Next, I go to the Catalog Tree Window and navigate to the location of the toolbox. In this situation, I can highlight the tool and choose the 'Description' tab and edit the metedata there. This, however, does not translate the edited metadata back to the same tool in the ArcToolbox window. If I open the tool in ArcToolbox window and highlight a parameter, I still get a description message of 'No description available'. So, it's like I have to edit the metadata twice. Once in the Arctoolbox window and then again in the Catalog Tree Window.

I've noticed this about the Description tab in ArcCatalog. It only honors selection in the Catalog Tree Window. If I a have tool selected in the ArcToolbox window, the description is still set to whatever is highlighted in the tree window, not the toolbox window.

Is this how it's supposed to work? I'm OK with only having to be able to edit metadata in the Catalog tree window, but the fact that metadata does not translate back to the same tool in ArcToolbox window is odd to me. Am I doing something wrong?


texts: Does Coordinate System of ArcMap data frame have to match at least one of its layers?

texts:
Does the Coordinate system of a data frame in ArcMap have to match the coordinate system of at least one of the layers in the table of contents?

I have found that if defined it will take on the coordinates of the previous layer.  If both undefined it will stay undefined?  I am not sure if the answer is yes or no.  I believe it is yes because I have found this excerpt.

  By default a data frame (the map display area) takes on the same
  coordinate system as that of the first spatial data layer added to it.
  It is important to realize this. But,assuming that all of the layers
  in your data frame have their coordinate systems defined properly, you
  can set the data frame to display your data in any desired coordinate
  system, regardless of the coordinate system of the layer you bring in
  first.

To me it is a bit confusing.  So if the coordinate systems are all defined then the answer would be no, I assume.  Can anyone clarify this for me?


texts: Merging and adding values based on ID column?

texts:
I have two shapefiles:

a grid containing just the grid "ID" column
a shapefile containing multiple entries from different road segments

The second shapefile represents emissions with a common grid "ID" as the one in the first shapefile (grids). This shapefile was a result of an intersection of the grid with the roads as seen below.

This is the attribute table of the first

This is the attribute table of the second shapefile

As seen, there are multiple entries for grid "ID" 229, 230, 231, etc. in the second shapefile.

What I want to achieve is a final grid shapefile such that each grid with a common grid "ID" has the total emissions for that grid.

How can I do this using QGIS?


texts: Is there a way to create an mxd from the command line?

texts:
Is there a way to run ArcMap in a batch script to create a basic mxd containing specified layers? Something like:

  "Path to ArcMap exe" -create "Path to new mxd" -include "layers to add"


texts: How to slow down mouse before modifyend with openlayers3

texts:
I have a map with a draggable marker (ol.style.icon in fact).
There is a strange issue:
If I try to move marker, it is like trying to reel in a marlin but after a zoomin/out, marker naturally slows down.
I tried to play with PixelTolerance but any difference.
Why and how to immediately slow down this marker, please?

var mark_style = new ol.style.Style({
    image: new ol.style.Icon({
    anchor: [.5, 48],
    anchorXUnits: 'fraction',
    anchorYUnits: 'pixels',
    opacity: 0.75,
    src: mark_path
 })
});

var iconFeature = new ol.Feature(position);
iconFeature.set('style', mark_style);

var vectorLayer = new ol.layer.Vector({
    style: function(feature) {
        return feature.get('style');
    },
    source: new ol.source.Vector({features: [iconFeature]})
})
var dragInteraction = new ol.interaction.Modify({
  features: new ol.Collection([iconFeature]),
  style: null,
  pixelTolerance: 10//?
});

dragInteraction.on('modifyend',function(f){
    var coordf = f.features.getArray()[0].getGeometry().getCoordinates();
    ol.View-view.setCenter(coordf);
    },iconFeature);

  map.addInteraction(dragInteraction);


texts: How is a "conformal latitude" used in ellipsoidal stereographic projections?

texts:
At my work, I've been doing some maintenance on internal software that claims to implement a Stereographic projection whose formulas were copied from an old system a long time ago.

From what I've seen, it seems to implement the "Oblique and Equatorial Stereographic" from EPSG (code 9809). But there is one step that I was not able to identify. The "R" value used in our code is calculated with the following parameters:

(phi_n, lambda_n) - null distortion coordinates
(PHI_nc, LAMBDA_nc) - conformal null distortion coordinates
(PHI_o, LAMBDA_o) - conformal projection origin coordinates
e - eccenttricity
r_eq - earth radius at equator

Given this, the formula to get R is:

a = 1 + cos(PHI_nc) * cos(PHI_o) * cos(LAMBDA_nc - LAMBA_o)
b = sin(PHI_nc) * sin(PHI_o)
c = 2 * cos(PHI_nc) * [1 - (sin(phi_n) * e/2) ^ 2] ^ 1/2
R = r_eq * cos(phi_n) * (a + b) / c

Does anyone know what is going on here? I've been searching and reading about stereographic projections, but I wasn't able to find anything like that. In fact, I've not seen any projections that use this "null distortion coordinate."


texts: R: How to write CRS information correctly to NetCDF?

texts:
As the title says, I'm struggling with writing the CRS of a multi layer raster stack to a NetCDF file using a custom package. This package uses the functions of the ncdf4-package. But the written .nc file doesn't have the CRS information which the source raster had. I've already read this thread and I tried to replicate it, tweaking the WriteNCDF function from above. I tried it with:
ncatt_put(nc, var.name, 'proj4', crs(data), 'double')

with ncas the to be written ncdf file and data as the rasterlayer.
This results in the following error:

Error in ncatt_put_inner(idobj$group_id, idobj$id, attname, attval, prec = prec,  :
ncatt_put: error, passed an attribute with a storage mode not handled.  Att name: proj4 Att value: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs Storage mode passed: S4 .  Handled types: integer double character

Afterwards I tried it with:
ncatt_put(nc, var.name, 'proj4', crs(data, asText = T), 'text')

This doesn't produce an error but when I open the written .nc file with raster::brick(), the CRS isn't found.
So how do I put CRS information correctly into a NetCDF file?
Minimal code to run:
library('ncdf4')
library('raster')
source('WriteNCDF.R')

r <- raster(crs="+proj=utm +zone=32 +datum=WGS84", xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- 1:ncell(r)
s <- stack(r, r/2)
x <- writeRaster(s, filename="test.nc")
crs(x)

z <- WriteNCDF(data.l = list(s$layer.1, s$layer.2),
               var.name = c('layer1', 'layer2'),
               var.unit = c('unit', 'unit'),
               file = 'test2.nc',
               utm = T)
zz <- brick('test2.nc', varname = 'layer1')
crs(zz)

The tweaked WriteNCDF.R can be found here.


texts: How to resolve "'NoneType' object is not callable" TypeError when calling rasterstats.zonal_stats from ArcPy script?

texts:
I'm attempting to workaround the apparent zonal statistics bug in ArcGIS 10.4 by using rasterstats.zonal_stats.

stats = zonal_stats(shapefile_full_path, rast_full_path,0,1,None,None,{'mean','min','max'},False,False,None,None,False,None,True )

In a standalone script, it works perfectly.

But when called from an ArcPy geoprocessing script by ArcGIS Server, I receive this error:

  File "C:\Python27\ArcGISx6410.4\lib\site-packages\rasterstats\main.py", line 28, in zonal_stats
    return list(gen_zonal_stats(*args, **kwargs))
  File "C:\Python27\ArcGISx6410.4\lib\site-packages\rasterstats\main.py", line 139, in gen_zonal_stats
    with Raster(raster, affine, nodata, band_num) as rast:
  File "C:\Python27\ArcGISx6410.4\lib\site-packages\rasterstats\io.py", line 237, in __init__
    self.src = rasterio.open(raster, 'r')
  File "C:\Python27\ArcGISx6410.4\lib\site-packages\rasterio\__init__.py", line 201, in open
    s.start()
  File "rasterio/_base.pyx", line 67, in rasterio._base.DatasetReader.start (rasterio/_base.c:2557)
  File "C:\Python27\ArcGISx6410.4\lib\site-packages\rasterio\vfs.py", line 30, in parse_path
    parts = urlparse(path)
TypeError: 'NoneType' object is not callable

It is a similar error to this question, but the answer doesn't apply (no RDBMS involved in my script):
ArcGIS Geoprocessing "None Type" Second time

There's no error in the from rasterstats import zonal_stats statement, and the error is raised in rasterio.

Based on this, it appears that this was successful at one point:
https://github.com/perrygeo/python-rasterstats/blob/master/docs/notebooks/ArcPy%20and%20rasterstats%20integration.ipynb

Any idea of how to resolve? Can it be resolved?

This is an ArcGIS 10.4.1 Server, but the python.executable is ArcSOC.exe (I suspect this has something to do with it):

Python executable:C:\Program Files\ArcGIS\Server\bin\ArcSOC.exe
Python version:2.7.10 (default, May 23 2015, 09:44:00) [MSC v.1500 64 bit (AMD64)]
Python hexversion:2070af0
os.__file__:C:\Python27\ArcGISx6410.4\Lib\os.pyc


texts: Add Interactive Flow Arrows to Story Map

texts:
Is it possible to add interactive arrows showing the movement of materials from one area to different sites?

The idea would be to create a story map showing the arrows (or something similar) at an overview and then the user would be able to click on each arrow separately which would direct them to a different page showing specific information about that facility.


texts: How to export attributes from more than one layer in QGIS3 to csv?

texts: I have two vector layers, one has geographical points and the other has polygons representing suburb boundaries. I want to export a list of all the points in a table with a column listing the suburb and postcode in which the points are located.


texts:
I am trying to download SRTM rester data in QGIS for a specific area. See image below:

I am using SRTM-Downloader plugin and when I click on download the following error pops up. Any idea about what the problem is?


texts: TileMill OGR Plugin: Failed to read GeoJSON data

texts:
I'm trying to load the following into a TileMill layer as the file build1.geojson

{
    "type": "Feature",
    "geometry": {
    "type": "Polygon",
    "coordinates": [[
        [752669,349453],[913955,-4919228],[4300965,-4865466],[4354727,-12499679],[-13010418,-12499679],[-13010418,-7123473]
        ]]
    },
    "style":{
        "fill":"red",
        "stroke-width":"3",
        "fill-opacity":0.6
    }
}

However, I'm getting the error shown in the title. I can't see anything wrong with the syntax so what is my error? Is it tilemill itself?


texts: QGIS 2.0 text delimited layer importing data as single row

texts:
I'm trying to import a text delimited layer using lat/lon, but the tool interprets the data as a single line.  I've tried saving my csv using several different programs (Stata, python editor, text editor, Excel), but nothing works.

Here is a sample of my data:

n,aqs_site_id,site_latitude,site_longitude
1,17-001-0006,39.93301,-91.40424
2,17-001-0007,39.91541,-91.33587
3,17-019-0004,40.12379,-88.22953
4,17-019-0006,40.12389,-88.24055

And my csvt file:

"Real", "String", "Real", "Real"

This is what the dialog box looks like when I try to upload.  You can see there is no data because everything is pulled into row 1:

Any help would be much appreciated!


texts: Is it possible to have a locked TOC in QGIS

texts:
I would like to customize QGIS for a custom use case where the user will work with a very simple customized UI and will just digitize features in existing layers.

Therefore the TOC should be locked, so that drag and drop reordering of the LayerTree will be blocked as well as inserting layer by drag and drop.

Is this possible?


texts: Calculating Service Area for each feature in feature class using ArcGIS Network Analyst?

texts:
I am trying to implement a python script that creates in a first step a Network Dataset and that in a second step gets a point feature class for which the tool should create for each feature in the feature class a service area. The first step is no problem, but the second...

My first idea was to use the cursors in python which should run through the whole attribute table. For each feature a temporary feature class should be created which should be used as input for the service area.

Is this the right direction or is there an more easy way to calculate for each feature individually a service area?


texts:
*Running QGIS 3.6, GRASS 7.6; have both Python 2 and 3

I'm trying to run a GRASS algorithm in QGIS, and none of the GRASS tools are working. GRASS also isn't working in the standalone version. I have tried both QGIS and QGIS with GRASS.

I tried changing the name of my Python 3.7 path as that has caused issues in the past, and got this error when running in QGIS:

  This algorithm cannot be run :-( The specified GRASS 7 folder
  "C:\Program Files\QGIS 3.6\bin\bin" does not contain a valid set of GRASS 7 modules.
  Please, go to the Processing settings dialog, and check that the GRASS
  7 folder is correctly configured

When I run GRASS independently (for 7.4 and 7.6), the command line opens with the following error:

File "C:\Python37\Lib\site.py", line 177
  file=sys.stderr)
SyntaxError: invalid syntax
Press any key to continue . . .

Running grass within QGIS (v.clean tool):

  The following layers were not correctly
  generated. [....] You can check the
  'Log Messages Panel' in QGIS main window to find more information
  about the execution of the algorithm.

There is nothing substantial in Log Messages, but looking through the GRASS log I found the same thing I was getting when trying to run GRASS standalone:

File "C:\Python37\Lib\site.py", line 177
file=sys.stderr)
^
SyntaxError: invalid syntax
Press any key to continue . . .

I checked my processing options: GRASS is activated and the docs are in the proper location.

I also tried looking for any errors around the 'syntax error' in Python37\Lib\site.py, as well as a couple of GRASS .bat files, but I don't know enough about using GRASS to actually tell if anything is wrong.


texts: Stacking label boxes in ArcMap label expression

texts:
Is there a way to automate getting a seperate text box for two different field names in the label expression in Arcmap?
Right now I have both fieldnames stacked in the same box. Without having to convert to annotation I would like to be able to have each field name in its own box.


texts: I'm trying to remove clouds from Landsat 5,7 and 8 data. I'm also trying to do this using only open source options, avoiding ArcGIS. I'm familiar with python and have some QGIS experience. Originally, I was trying to use the Semi-Automatic Classification model in QGIS, but it is giving me a lot of problems, and doesn't always output actual classifications (it often just returns a uniformly valued square...). Can anyone point me in the direction of a better option?


texts: Error during MGRS conversion for large latitudes

texts:
I'm using Nasa's Worldwind java API to convert MGRS strings to lat/lon. However, when my latitude is greater than 85 degrees or less than -85 degrees I get an error.
 java.lang.IllegalArgumentException: Error during MGRS conversion

Here is my input

MGRSCoord.fromString("ZAD 00000 66856", null);

I'm creating the MGRS from worldwind was well.

Double latitude = 87.0;
Double longitude = 0.0;
Angle latitudeAngle = Angle.fromDegreesLatitude(latitude);
Angle longitudeAngle = Angle.fromDegreesLongitude(longitude);
String coordString = MGRSCoord.fromLatLon(latitudeAngle,longitudeAngle).toString();
MGRSCoord.fromString(coordString, null);

My guest guess is this is frdefault Globe into the conversion by using null, but I tried using actual Globes, but it didn't work, and frankly I didn't really know what to set it to.


texts: Copy and paste from Microsoft Excel to ArcGIS 10.1 is not showing pivot table?

texts: When I work on my pivot table in MS Excel 2013 and made a copy of it and paste to my ArcGIS 10.1. Why does it not showing ? Is there a work around to make it work ?


texts: Calculate the mean and the use of list lead to lost time property

texts:
I want to output charts of mean NDVI per month, but calculate the mean and the use of list lead to lost time property. And I don't kown how to copy property.

//渲染Layer显示颜色
var vis = {min: -0.2, max: 1, palette: [
  'FFFFFF', 'CE7E45', 'FCD163', '66A000', '207401',
  '056201', '004C00', '023B01', '012E01', '011301'
]};

//定义图表中不同作物类型的颜色
var COLOR = {
  LINDI: 'ff0000',
  GENGDI: '0000ff'
};
//定义感兴趣区
var lindi = ee.Feature(
    ee.Geometry.Rectangle(115.105882, 35.430421,115.106077, 35.430255),
    {label: 'LINDI'});
var gengdi = ee.Feature(
    ee.Geometry.Rectangle(115.104388, 35.430467, 115.10501, 35.429855),
    {label: 'GENGDI'});

//将作物感兴趣区合并
var cropRegions = new ee.FeatureCollection([lindi,gengdi]);

//去云
var s2_rmcloud = function(image) {
 var quality = image.select("QA60").unmask();
 return image.updateMask(quality.eq(0));
};
//compute ndvi
function NDVI(image) {
 return image.select().addBands(image.normalizedDifference(['B8', 'B4']).select([0], ['NDVI'])
 .copyProperties(image, ['system:time_start']));}
var start_list=new Array('2018-1-1','2018-2-1','2018-3-1','2018-4-1','2018-5-1','2018-6-1',
      '2018-7-1','2018-8-1','2018-9-1','2018-10-1','2018-11-1','2018-12-1');
var end_list=new Array('2018-2-1','2018-3-1','2018-4-1','2018-5-1','2018-6-1','2018-7-1',
    '2018-8-1','2018-9-1','2018-10-1','2018-11-1','2018-12-1','2019-1-1');

var test=[];
for (var i=0; i<start_list.length; i++)
{
  var s2= ee.ImageCollection('COPERNICUS/S2')
 .filterDate('2018-1-1', '2019-1-1')
 .filterDate(start_list[i], end_list[i])
 .filterBounds(cropRegions);

var s2 = s2.map(s2_rmcloud)

//筛选s2数据

for (var i=0; i<start_list.length; i++)
{
  var s2= ee.ImageCollection('COPERNICUS/S2')
 .filterDate(start_list[i], end_list[i])
 .filterBounds(cropRegions);

var s2 = s2.map(s2_rmcloud)
//计算每副影像的NDVI并制作数据集
var ndvi = s2.map(NDVI);
//Map.addLayer(ndvi, vis, 'NDVI_'+i);

test.push(ndvi.mean().copyProperties(ndvi));
//Map.addLayer(test,vis,'NDVI_average');
//print("test", test)
}

  var test = ee.ImageCollection(test)
var ndviTimeSeries = ui.Chart.image.seriesByRegion({
  imageCollection: test,
  regions: cropRegions,
  reducer: ee.Reducer.mean(),
  band: 'NDVI',
  scale: 10,
  xProperty: 'system:time_start',
  seriesProperty: 'label'
});
ndviTimeSeries.setChartType('ScatterChart');
ndviTimeSeries.setOptions({
  title: 'Sentinel-2数据作物时间序列NDVI变化',
  vAxis: {
    title: 'NDVI'
  },
  lineWidth: 1,
  pointSize: 4,
  series: {
    0: {color: COLOR.LINDI},
    1: {color: COLOR.GENGDI},
  }
});

Map.addLayer(cropRegions, {color: COLOR.LINDI},'ROI');
Map.centerObject(cropRegions, 8)
// Map.addLayer(lindi, {color: COLOR.LINDI});
// Map.addLayer(gengdi, {color: COLOR.GENDDI});
print(ndviTimeSeries)


texts: Doing a viewshed analysis from INSIDE a building?

texts:
I'm currently trying to do a viewshed (spatial analyst), using ArcGIS.
The datasets being used are observation/observer points located WITHIN buildings, and a DSM of the surrounding area and the buildings.
To force the observation points inside the buildings, i have set the SPOT attribute to values lower than the buildings themselves (i.e. the observer is placed on the floor inside the building, instead of the roof).

Running this analysis is not possible, since it gives me the following error:

ERROR 010067: Error in executing grid expression. No valid Observation points present.

Reading through an answer from another thread on this website, the error occurs because the SPOT value is lower than the input DSM value in the same location (Error in viewshed (arcgis)).

I was wondering if there are any smart solutions/workarounds that mean i dont have to place the observers on top of the buildings?


texts:
I have an .evf file and it is not opening in QGIS.

How can we reproject the vector layer (specifically if we want to convert it from lat/lon to UTM projection)?


texts: Debugging syntax error (invalid syntax) from arcpy.SelectLayerByAttribute_management()?

texts:
I'm writing a script that will create a list of values from the field "Year_" using search cursor, and then a set() of unique values from the list. I then want to iterate through the set making a main folder with the current iteration value for "Year_" followed by the creation of a selection variable also based on the value of the current iteration. From this selection I want a second list made based on the values in the "Customer" field and then a set of unique values from the list. Finally, the script should make sub folders of unique customer names based on the unique values in the set:

    # Import modules

    import arcpy, os , sys, string

    # Create environmental variables

    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = r"C:\data.gdb"

    # Set variables

    input = r"C:\data.gdb\layer"
    yearField = "Year_"
    custField = "Customer_Name"
    rootpath1 = r"S:\Shared\Orders"

    # Use Search Cursor to make list of years and then make set of years with unique values
     yearList = [row[0] for row in arcpy.da.SearchCursor(input,(yearField))]
     yearSet = set(yearList)

     #Make main year directories
     for year in yearSet:
         rootpath2 = os.mkdir(os.path.join(rootpath1,year) #establishes root  path for subfolders while making main directories
         selection = arcpy.SelectLayerByAttribute_management(input,year)#creates variable based on current selected year in input layer
        for customer in selection:
            customerList = [row[0] for row in arcpy.da.SearchCursor(input,(custField))] #Use Search Cursor to make list of customers and then make set of   customers with unique values
            customerSet = set(customerList)
            for name in customerSet:
                os.mkdir(os.path.join(rootpath2,name) #makes unique customer sub folders

So if the values for years are 2010, 2010, 2011, 2012, 2013, 2013 in the field "Year_" a list will be made and the resulting set would be [2010,2011,2012,2013].
The first iteration would find 2010 and make the pathway r"S:\Shared\Orders\2010".

While in this iteration, 2010 would be used to select by attributes. If customers Joe, Joe, Ann, Bob, Mary, Mary, Mary were found a list would be made and then the set [Joe, Ann, Bob, Mary]. Finally a for loop will iterate through the set and make the pathways:

r"S:\Shared\Orders\2010\Joe"
r"S:\Shared\Orders\2010\Ann"
r"S:\Shared\Orders\2010\Bob"
r"S:\Shared\Orders\2010\Mary"
However, I get a an error "Failed to run script - syntax error - invalid syntax at the line where I make the selection variable.


texts: ArcGIS Online Routing - Multiple route end stops?

texts:
Is it possible within the ArcGIS Online analysis - Plan Routes, to select multiple route end stops?
By default you can can select a stop to start the routing (depot), the stops to visit, and an end point. This is fine if for example I have a school bus leaving a depot, picking up 50 pupils, and arriving at a single school.
What I am trying to do is arrive at multiple schools. So it will start at a depot, pick up 50 pupils, and drop off at say 3 different schools.
I can't seem to see any options within AGOL for this, does anyone know if it is possible?


texts: how to create a shapefile for 2 million points with geotools

texts: I have to create a shapefile with the geotools java library for approx.  2 000 000 points and whatever strategy I apply I seem to run into out of memory exceptions as well as everything running super slow. Anybody done that before and got it to work?


texts: How do I prevent Arcmap from verifying SQL expressions when using a plugin dataset?

texts:
I've written a plugin workspace for ArcGIS.

I'm able to compose a where clause containing syntax that my implementation of IPluginDatasetHelper understands, assign it to IQueryFilter.Whereclause, and get the expected results from ITable.Search.

However, if I copy and paste the same whereclause into into the layer query Definition expression editor, it fails to verify.  Arcmap complains that it is not SQL.  In this case, the expression never makes it to my plugin.

Is there someway I can escape my definition expression in Arcmap so that it lets the plugin validate he expression?


texts: Is possible to use ArcView select and query data from Geoserver WFS?

texts:
I am looking at implementing a GIS solution using GeoServer.

If I implement a GeoServer solution is possible using ArcView via WFS to select and query data?  From the research I have done it looks like there is a client plugin cost of around $1500 (Data Interoperability extension).

As most of our clients use ArcView it is not practical to ask each of them to purchase this extension.


texts: QGIS: where did these user-defined coordinate systems come from?

texts: At 3.16.1, when I view the Project Properties > CRS page, I notice 25 user-defined coordinate systems that I don't ever remember creating (within the screenshot red box, below).  I can only assume that they arrived by default when I installed QGIS.  Why are they there?  Can I get rid of them?


texts: Calculate river length in each subcatchment

texts:
Using ArcGis 10.4

I have suceeded in calculating "maximum upstream river distance to its reach per subcatchment" per subcatchment with the help of this forum (see Calculate river distance to headwater for every pixel?).

Now I additionally need the length of the river stretch within each subcatchment. My idea was to just calculate the minimum upstream distance to its reach and then deduct this from the previous maximum distance. But it doesnt work, I run into multiple problems when joining them later. So I was thinking, zonal statistics as table should do the trick. But it doesnt work either.

Any other idea how I could calculate that?


texts: Reading feature class in file geodatabase using R?

texts: I have a feature contained in a geodatabase that is larger than 2GB as an exported shapefile.  I need to run an extract function in R to attribute the polygons with data from a raster file.  Exporting the feature as a table is not a solution. How can I read feature classes contained within an Esri file geodatabase?


texts: How do I get unique names for model builder outputs?

texts:
I have set up a model that erases designated restricted areas such as Flood zones, SSSI, national parks etc, each erase is linked from the output to the next erase going step by step through all the different restrictions, 12 restrictions altogether.

My first input is called Usable_Parcels_10m_Rail_Line_Editing, what I would like is each erase output to be named as Parcels_10m_Rail_Restriction with the restriction part replaced by the specific restriction name from the restriction feature class named as SSSI_Merge, Nat_Parks_merge etc excluding the merge bit?

Is this possible?


texts:
I'm generating an atlas in QGIS and adjusting it to only show the legend when the items are inside the current atlas feature. The problem is that when the items are not inside the atlas features it still appears as a small white box that corresponds to the background of the legend. I don't want to set the background to transparent as I have a satellite image underneath and it would be difficult to read the legend.

Is it possible to make the legend box completely disappear when my items are not in the atlas feature?


texts: How to close all (or distinct) attribute table with PyQGIS

texts:
I have a QGIS (2.18.0) plugin which updates a join between to layer (spatial and non-spatial layer).

QGIS crashes when I run the plugin and the attribute table of the joined layer is open. There is no crash when the attribute table is closed or an attribute table of another layer is open.

The last two commands that the plugin runs are:

....
PUNKTE_OS.addJoin(joinObject)

PUNKTE_OS.triggerRepaint()
....

If there is no other solution, how can I close the attribute table of the joined layer before the new join is set?

I know, there is the showAttributeTable() method. Is there also a "closeAttributeTable()" method?

UPDATE:

It works when I close the specific attribute table in the def __init__(self, parent): function on Plugin startup. It doesn't work when I close the specific attribute table in the function where I set the joins (setJoin(self)).

Here a sample of my code:

class mask_n_functions(QDialog, Ui_Dialog):

def __init__(self, parent):
    QDialog.__init__(self, parent)
    self.setupUi(self)
    ...
    attrTables = [d for d in QApplication.instance().allWidgets() if d.objectName() == u'QgsAttributeTableDialog' or d.objectName() == u'AttributeTable' ]
    for x in attrTables:
        if 'attr_table' in x.windowTitle():
            x.close()

def setJoin(self):
    prjpath = QgsProject.instance().fileName()
    ...


texts:
I am trying to get GDAL running under the following environment:

Win 10
Miniconda3 (this is just Anaconda with fewer modules loaded by default)
Python 3.6.6

I installed GDAL with
conda install -c conda-forge gdal

When I try to import osgeo.gdal, I get this error:

(base) C:\Users\askov>conda activate geoprocessing

(geoprocessing) C:\Users\askov>python
Python 3.6.6 |Anaconda, Inc.| (default, Jun 28 2018, 11:27:44) [MSC v.1900 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> from osgeo import gdal
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\askov\AppData\Local\conda\conda\envs\geoprocessing\lib\site-packages\osgeo\__init__.py", line 21, in <module>
    _gdal = swig_import_helper()
  File "C:\Users\askov\AppData\Local\conda\conda\envs\geoprocessing\lib\site-packages\osgeo\__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)
  File "C:\Users\askov\AppData\Local\conda\conda\envs\geoprocessing\lib\imp.py", line 243, in load_module
    return load_dynamic(name, filename, file)
  File "C:\Users\askov\AppData\Local\conda\conda\envs\geoprocessing\lib\imp.py", line 343, in load_dynamic
    return _load(spec)
ImportError: DLL load failed: The specified module could not be found.

I saw some articles (for example) suggesting I install GDAL from GIS Internals, which I did. The GIS Internals version I used was gdal-203-1900-x64-core.msi. I also updated the PATH, GDAL_DATA, and GDAL_DRIVER_PATH, but that had no effect. I checked those variables in my Conda environment, and they are different from the system settings. Even after I manually set them (via the DOS set command), I get the same result.

Another article suggested installing the krb5 module, which I installed via conda, with no change in the result.

How can I get gdal working in from a clean python environment (i.e.: not via QGIS or OSGeo4W)?


texts:
I am working on a provincial annual health bulletin.

I have a .shp file that has GIS data about the smallest divisions called "GN Area". Two or three GN Areas form one MOH area. I want to get the MOH area details displayed in the map. How can I do it? I am using QGIS.

Is there any way where GN Area data exported to excel and add the relevant MOH area and import it back to Quantum GIS and display MOH Areas rather than GN ares?

If that is not possible with Quantum GIS, is there any other method to do it?


texts:
I am trying to georeference a pdf into my QGIS project.

I am doing the exact same way as I have done many many times (same transformation settings, number of reference points...). After doing the georeferencing, the raster shows up in my layer panel but it is not visible. Whether I hide all other layers, change transparency, zoom to layer...the layer is there but invisible.
I'm using 3.4.10 (Madeira) which is the version I have been using for months

What am I doing wrong?


texts: What are the corner coordinates of a global map in meters?

texts:
We all know that the corner coordinates of a global map in degrees are -180, 90; -180, -90; 180, 90; and 180, -90 (top left; bottom left; top right; bottom right).

So what are the corresponding coordinates in meters in arcgis using D_WGS_1984 datum?

I need to know these coordinates to define projection for my MODIS blue marble image.


texts: How to update a hosted feature service on ArcGIS Online using the Python API for ArcGIS Pro?

texts:
I'm trying to follow the steps in the Esri blog post Updating your hosted feature services with ArcGIS Pro and the ArcGIS API for Python but am running into an error.

To test this I published a simple stand-alone file GDB into an ArcGIS Online web map using ArcGIS Pro's UI tools. This creates a Service Definition, WebMap and associated Feature Layer on ArcGIS Online.

I then made a simple symbology change in ArcGIS Pro and attempted to push the updates to ArcGIS Online using the script. The script runs correctly as far as these lines:

sdItem = gis.content.search(query="title:"+ sd_fs_name + " AND owner: " + user, item_type="Service Definition")[0]
sdItem.update(data=sd)

I can confirm from the timestamps on ArcGIS Online that the Service Definition has been updated, but the script is falling over when trying to update the associated Feature Layer. The script contains this line:

fs = sdItem.publish(overwrite=True)

This gives me the error message

  RuntimeError: User cant overwrite this service, using this data, as
  this data is already referring to another service

I can confirm that the data doesn't actually refer to another service since it's coming from a stand-alone file GDB that isn't used anywhere else. I'm using the same ArcGIS Online administrator's login within ArcGIS Pro and the Python script so I'm definitely the owner of the layer.

How can I overwrite the feature layer using the latest changes in the ArcGIS Pro document?


texts: Old Windows Installers for python library fiona

texts:
Due the the recent bug as presented here for fiona 1.1, i am obliged to use a version of fiona in 1.0.x on windows 8.1. however, i can't find older versions of the installer at here.

Are there any sources for the archive installers, specifically for the 32bit fiona (1.0.x) ? thanks in advance!

P.S.

I think it will apply to ALL other previous libraries, and not just fiona


texts: What is the meaning of this linear scale of angles on this 1951 map of Iceland from the U. S. Army Map Service (AMS C762)?

texts: On the U.S. Army Map Service (AMS C762) 1:50,000 maps made of Iceland 1948–1951, there is a scale marking at the top which shows degrees, I suspect angles.  In the image below it shows 25°–29°.  It increases to the left, which suggests it might refer to degrees longitude west, but the numbers are not correct; the map covers 18°W–18°30'W.  It's not a magnetic declination either, which is shown elsewhere on the collar in the usual way.  Other maps in the series have a similar scale in the same place, but the numbers may be different.  What does this angle scale mean?


texts: I need to lock a file geodatabase (or only a single Mosaic dataset) at the beginning of a python script (and release it in the end).. I looked and looked and couldn't find any explicit arcpy command to do this.. Is it even possible to tell the gdb to lock itself? Or is the locking mechanism works autonomously?


texts: How to average # of parks per child within 1 sq mile and average parks per child in whole county?

texts: I have a shapefile of census tracts within Cook County, IL that has been combined with a shapefile of every park in the cook county area. It displays as a graduated color scheme of number of parks per tract. On top of that, I have a dot layer from the that displays the number of children per census tract in the county. Each dot represents 500 known children. I'm trying to form a poster that displays the relative access that children have to parks in the county. As the project stipulates having supporting graphs, tables and statistics, I would like to display the average amount of parks available to a child in contrast to the average amount of parks per child in the county. I will attach a snippet of the map for reference.


texts: Is there a way to manually separate the distance between 2 lines so that it is represented better?

texts: I work with powerlines and there are too many lines quite close to each other making it look line one thick line when zoomed out . Obviously they do separate out when zoomed in but most of the time we need to work at a zoomed out level


texts: qgis2compat can't be installed in Yakkety

texts:
I want to use QField and Qfield sync. QField depends of qgis2compat and I can't install it.

I have QGIS 2.14.5 and Ubuntu 16.10 Yakkety and can't install qgis2compat plugin because the lack of QtWebKit because python-qt4-webkit has been removed.

Is there a workaround?


texts: Clipping raster data with shapefile in ArcPy

texts:
Am trying to clip a raster data to the extent of the input shapefile but I keep getting an the following error:

Traceback (most recent call last):
File "C:\Python27\dealingwithraster.py", line 40, in <module>
arcpy.Clip_management("outputRaster","#","outputclip","input_studyArea","ClippingGeometry","MAINTAIN_EXTENT")
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 13771, in Clip
raise e
ExecuteError: ERROR 000622: Failed to execute (Clip). Parameters are not valid.
ERROR 000800: The value is not a member of ClippingGeometry | NONE.
Failed to execute (dealingWithRaster)

Below is the my code:

import arcpy
input_raster= arcpy.GetParameterAsText(0)
input_studyArea= arcpy.GetParameterAsText(1)
outputRaster = Con(inRas1>=0 ,inRas1,1)
outputRaster.save(outputgdb + rasterName)
arcpy.Clip_management("outputRaster","#","outputclip","input_studyArea","ClippingGeometry","MAINTAIN_EXTENT")


texts: Create a Raster file out of a ASCII grid file?

texts:
I'm using a special system to do geomagnetic surveys. The software let's me output the information as image files, ASCII grid or in its own proprietary file format.

Now I want to get its information in my GIS (using QGIS). I can of course import the exported image files but I would love to import the raw data to play around with the visualisation a bit more, not depending on the original software.

The files look like this:

X   Y   Messwert
995.00  1000.00 0.00
995.00  1000.05 0.14
995.00  1000.10 0.28
995.00  1000.15 -0.07
995.00  1000.20 -0.42
995.00  1000.25 -1.26
995.00  1000.30 -2.17

and so on..... (looooong file ;) )

Every 0.05 cm movement on X and Y contains a data value, positive and negative.

I want to convert this data into a raster file, a pixel for every data value.
My goal is to tune the visualisation comparable to a DEM TIFF file in my QGIS project, without having the problem to export this through the original software every time.

What would be the best way to do this?

I think GDAL is the program to use I need some help. Perhaps there is even a way to do this in QGIS?

Update:

So I finally imported a fraction of my data, sorting everything to Y.
Saving everything as TIFF wasn't a problem as well. Now my next step is to get this spatially correct data (in terms of length) into my project. The coordinates in the file are just a local project oriented coordinate system.

Georeferencing the created TIFF wasn't a big problem but it results in a little annoying problem. after georeferencing my perfect square gets rotated a bit, resulting in big nodata areas.

My data also contains positive as well as negative data and even zero is important.

I couldn’t find a way to get this nodata area to disappear, QGIS  georeferencing gives it a value that is contained in the data areas as well.
if I set this to transparency my raster files gets some annoying holes.


texts: QGIS Table Manager plugin does not appear

texts: I installed Table Manager Plugin, but cannot find it under the plugin menu. Anybody know where to find it?


texts: Clipping and exporting sections of raster from multi-feature shapefile using ArcMap

texts: I have a multi feature shapefile and a slope raster and I'm trying to export images (tif, png whatever) for each feature showing the feature boundary and the slope raster contained within it with a field from the shapefile as the filename. I am dealing with a lot of features (>100,000) so efficiency is key, I tried data driven pages and it failed every time.


texts:
I'm developing a plugin in Python for QGIS.

I would like the user to select a particular area in a list (let's say a country list) and zoom in the map on the corresponding area.

Actually I'm doing it 'ugly':I create a layer (from a WKT), zoom to this particular layer, and delete this layer...

Is there a way to do it properly?


texts: Creating surface layer from Points in ArcScene

texts:
I work in ArcScene where i have some buildings in 3D, on the facades of these buildings in have points of solar intensity. I would like to create a continuous surface of solar intensity for each facade. I tried Interpolation with kriging but doesn't seem to work,

Does anyone have an idea how can i achieve this?

I join a picture for a better understanding.


texts:
I am working on creating a change detection script for QGIS3, that detects the differences between two DEMs and highlights the changes. Everything is working well but I have been unable to get the 3D view to work using pyqgis, I have it working using the GUI.

I want to display my output rasters in 3D view using PyQGIS, like this:

I have looked into the 3D Library and the Python API but I have not been able to get any of this to work.

I tried using the Qgs3DmapScene Class but the QGIS python console does not recognize it.

How can I use the 3D View to display layers using PyQGIS? is there any example or sample of this somewhere?

I am using QGIS 3.4.4 on Windows with Qt 5.11.2


texts: Filtering by folder or authoritative maps when using gis.content.search

texts:
I use the code in Python
for item in gis.content.search(query="owner:MyUsername", item_type='Feature Service')
to go over the feature classes uploaded on my organization's AGOL page. I want to filter out layers that are

Not View Layer
Authoritative
Are in a specific folder

Yet I had no success. I tried things such as
folder='MyFolder, status=authoritative
and I fail. Any solutions or workarounds on how to filter for these are welcome.


texts: Creating signature file for classification in ArcGIS Desktop?

texts:
I am trying to classify a raster image with the Maximum Likelihood Classification tool.

I understand from the ArcGIS Resource Center I need to input a signature file in the second field to tell the tool how to classify it. But I don't see how to make the signature file (see figure below). I think the Input raster bands field just needs the raster I am classifying.

What goes into the Input raster or feature sample data field?

The raster I have entered in the Input raster bands field is the same file I want to classify.


texts:
I need to print a map (with Arcgis 10.0) that has an extension of 1200 kmq that formally has a scale of 1:25000.
Simply I need to print it with that specified scale (1:25000)  but do not know how.

I've tried to looking on online guides, arcgis guides etc, and tried many options on print setup without success.

Some suggestions? Thanks!


texts: Keep selection order when Copy & Paste From Attribute Table in Qgis

texts: When copying some selected points listed in the Attribute Table from QGis to Excel the order seems to get lost. Is there a way to keep the selected points in the order they were originally selected?


texts: How to install multiple versions of PostgreSQL in same os Windows 7 64-bit?

texts:
I have already installed PostGIS along with OpenGeo suite and it is the latest version. Now I want to run pgrouting but the latest stable version is not compatible to the latest PostgreSQL versions.

Can I downgrade the PostgreSQL from 9.2 to 8.3?

I can't try some unstable versions of pgrouting that are compatible with PostgreSQL 9.1/9.2 because, there are issues in shooting star algorithm... or can I install another version parallel to this?

I tried uninstalling OpenGeo suite completely and installed only PostgreSQL 8.3 and PostGIS 1.4 but there are more and more errors I am getting and I am afraid to solve those as it's for my mini project and I have only 3 days to prepare documentation...

So someone please help! I am trying but I am confused with these kind of installations.


texts:
Given a vector of two points with direction from point 1 to point 2:
latitude_1 longitud_1 latitude_2 longitud_2
41.378118   2.169557   41.389474    2.156421

Find points in the dataset that are in the same "direction":
Latitude Longitud
41.383093   2.181116
41.373258   2.159358
41.385252   2.168779
41.390692   2.148911

In the image below, the 3 red points in front of the vector would be in the same direction.
For direction we can divide 360 degrees in 6 equal parts. A 4 part split is also good.
Expected output would be:
latitude_1 longitud_1 latitude_2 longitud_2   Number of points
41.378118   2.169557   41.389474    2.156421  1

Things I have tried with a 4 quadrant split:
def relative_location(row):
    if row["lat_dif"] > 0 and row["long_dif"] > 0:
        return "North East"
    elif row["lat_dif"] < 0 and row["long_dif"] > 0:
        return "South East"
    elif row["lat_dif"] > 0 and row["long_dif"] < 0:
        return "North West"
    elif row["lat_dif"] < 0 and row["long_dif"] < 0:
        return "South West"
    elif row["lat_dif"] == 0 and row["long_dif"] < 0:
        return "West"
    elif row["lat_dif"] == 0 and row["long_dif"] > 0:
        return "East"
    elif row["lat_dif"] > 0 and row["long_dif"] == 0:
        return "North"
    elif row["lat_dif"] < 0 and row["long_dif"] == 0:
        return "South"
    else: return "Same point"

df["relative_location"] = df.apply(relative_location, axis=1)


texts: Why font size doesn't set for text in MapBasic?

texts:
I have annotations layer with that was imported from ArcGIS and want to set font in MapInfo.

ArcGIS annotations contains its settings inside annotations objects (in Shape field, as I guess) but this information duplicates in usual, non-gis string and integer fields (FontName, FontSize, Bold etc).
I'm read font options from that fields and try to set text font:

Dim fs As Integer
Dim altering_obj As Object
' ...
Fetch Next From annotations
fs = annotations.bold + 2 * annotations.italic + 4 * annotations.underline
altering_obj = annotations.obj
Alter Object altering_obj
    Info OBJ_INFO_TEXTFONT, MakeFont(annotations.fontname, fs, annotations.fontsize, BLACK, -1)

Font name, style (bold, italic etc), color are applies but font size doesn't.
May be I'm do wrong something? How can I change font size?

MapInfo versions is 10.5.2, and also I tried in very old 6.5.
Both versions has same behavior.


texts: How to display different SVG icon in the same layer in QGIS

texts:
I am working with a ponctual layer in QGIS 2.12.3 (Lyon). I have a determine a few points in my layer. Each point represents a different activity (ex : climbing, hiking, etc). In order to have a proper symbology, i have imported a few SVG icon.

At this time, i want each of my point to be represented by the corresponding activity. I am unable to do so. Could you help me figure it out?

Here is an example of what i am looking for

Each point with pointed with an arrow = Hiking
Each point with a star = Rock climbing
Each point = Kayaking


texts: GDAL - Perform Simple Least Cost Path Analysis

texts:
I am investigating methods to perform a simple least cost path analysis with gdal.  By simple, I mean using the slope of a dem as the only cost factor.

I would prefer to do using the python or .net bindings, but will take anything.  Can anyone suggest any good tutorials or the like?


texts: How do you set the position of a plugin window when you start the plugin?

texts:
I'd like to set the position of a plugin's window or dialog relative to the map canvas frame.
For example, if I start a plugin, I want it to appear somewhere in the left upper corner of the map canvas.

I can use absolute screen coordinates (choosing appropriate values for the first two parameters of the setGeometry method of a dialog), but that's not a flexible solution as users might adapt the size of the map canvas frame according to their wishes.

FYI: I use the Plugin Builder plugin to take off...

Based on @Matthias Kuhn's concise answer, I found a solution which I share here below. Note that the code corresponds to the script that is produced by the Plugin Builder plugin. Also note that the argument of mapToGlobal() must be a QPoint.

def run():
    self.mapcanvas_upperleft = self.iface.mapCanvas().mapToGlobal(QPoint(0, 0))
    offset_x = 24
    offset_y = 43
    pos_x = self.mapcanvas_upperleft.x() + offset_x
    pos_y = self.mapcanvas_upperleft.y() + offset_y
    self.dlg.setGeometry(pos_x, pos_y, 0, 0)

    # show the dialog
    self.dlg.show()
    # Do this a first step... then delete
    print self.dlg.pos().x()
    print self.dlg.pos().x()

In order to define the offsets, you can print the upper left corner of your dialog before you change its position with setGeometry. The pos() method returns a QPoint.

(hint: in QGIS, click Plugins > Python Console to see the print output)

As such, you know the difference between the global coordinates of 1) the upperleft corner of the map canvas and 2) the upperleft corner of the plugin dialog, which then helps you to define the offset you prefer.


texts: I used 1.7.4 qgis release on Windows 7, by which it was possible to manage ECWfiles; after passing to 1.8 release, I couldn't deal with an ECW file: it tells me that .ecw is not a supported kind of file on QuantumGis. What can I do? Is there anyone that could help me?


texts: Aggregate functions not allowed in trigger

texts:
I need help with a postgres/postgis trigger.

I have a table of retaining walls, each of which is associated with the adjacent road, except where the wall abuts a junction, where it will be split and associated with the adjacent road. However it retains the same ID & name.  Each wall polygon has a width of 1m and so its length equals its area.  I have a field in the table which records the area/length and a working trigger to populate it.
What I need is for another trigger to sum the area fields of all of the relevant records and write it to a different table.

The code from the trigger that produces the aggregate functions error is

g:=NEW.gid;
  UPDATE roads.rds_retaining_wall_register a
SET length = sum(b.poly_area)
FROM gazetteer.lsg_special_designations b
WHERE a.refno = b.refno
AND b.gid=g;

After searching for the answer, I think I need to use a select statement to do the sum and then use that to do the update, but being new to postgres, I'm stuck with the syntax of how to do it.


texts: Seeking publicly available road shapefiles with elevation fields

texts:
I need to trace paths over a road network in several places in the Mid-Atlantic region of the United States. I am wondering if there are any publicly available detailed road datasets with elevation fields for modeling nodes where no turns are possible (such as bridges/overpasses). I know that there is the Esri North American Detailed Streets dataset with its F_ZLEV and T_ZLEV fields, but that is now deprecated.
Are there any more recent equivalents of that dataset?
Do Census TIGER files contain the data for this?
Will I simply have to look for this on the GIS sites of every state in question?


texts: Accessing "valueField" property of layer's symbology using ArcPy?

texts:
This has been really hard to debug, as I had no idea where to look at first, but I have made an attempt to access the "valueField" property of a layer file.

This layer file is of the UNIQUE_VALUES symbology type (although I knew it, also successfully tested with code), and thus must have a read/write valueField property associated with it according to the Help (http://resources.arcgis.com/en/help/main/10.2/index.html#//00s30000005s000000).

However, as soon as I try to access it like:

valField = lyr.symbology.valueField

the script stops dead in its track. Not knowing exactly how to access it (there is only a write example in the Help), I also tried:

valField = lyr.symbology.valueField()

to no avail.

I then noticed in PyScripter, after further configuration of Pyscripter to allow auto-completion with arcpy, that there is no valueField property / method shown... See the first screenshot (no, it's not down the scrollable list). The second screenshot shows I have configured PyScripter properly, as you can see the arcpy auto-completion at work.

What am I missing? I am able to successfully access other properties like labelClasses.

ArcGIS 10.2.0 by the way...


texts: Is it possible to grow QGIS symbol size over time?

texts:
I've got a data set of people movements within my municipality. (between points of places)

I've managed it with the plugin for QGIS 'time manager'. Over time people come and go, and based on a count of movements of a specific point at a specific time (time manager), i want to let my point symbol grow.

Is this possible with QGIS rule based symbol?


texts: In Arcpy, Erase function error: "XML being loaded could not be parsed"

texts:
A error in my previously-smooth-running arcpy script has me confounded. The error occurs where the script selects individual features in two different feature classes, and then performs an Erase function.

It seems to be working fine for the first one or two files that it erases, but then
stops with "Error 000210: Cannot create output."

Using arcpy.GetMessages gives greater detail...

  ERROR 000210: Cannot create output
  The XML being loaded could not be parsed.
  [GDB_Items] The XML being loaded could not
  be parsed. [XML document must have a top level element.
  XML document must have a top
  level element.

The only useful info online I've found suggested restarting the script each time it fails, which works - in that I restart, and the script moves on to the next item in line to be processed .

This is (most of) the offending script, and I've highlight where the Error is coming from:

import arcpy
from arcpy import env
import time
from time import strftime

#Set overwrite output
arcpy.env.overwriteOutput = True
print arcpy.env.overwriteOutput

donutlog = "C:\\BOEM\\Donuts\\donut_log.txt"
workspace = "C:\\BOEM\\pt_buff2\\"
fileList = []
buffersd4 = []
buffersd3 = []

buffs = workspace+"buffer.gdb"
directorypath = buffs

#Insert test to see if ID was done around here:
finishedIDs = {}

FClist = []
with open(donutlog) as f:
    for line in f:
        if line.startswith("Executing"):
            FCname=line.split(' ')
            FClist.append(FCname[2][:-4])

for dirpath, dirnames, fgdb in arcpy.da.Walk(directorypath,
                                              datatype="FeatureClass",
                                              type="Polygon"):
    for fc in fgdb:
        fc_split=fc.split('_')
        StartTime=time.clock()
        if fc.endswith("d4"):
            if fc in FClist:
                print fc + " has already been looked at!"
            else:

                #list features by APR16_ID in that fc
                APR16_ID4 = [row[0] for row in arcpy.da.SearchCursor(buffs+"\\"+fc,['APR16_ID'])]

                # make feat layer of fc_endswith_d4
                in_feat4 = buffs+"\\"+fc
                out_layer4 = fc+"_lyr"
                arcpy.MakeFeatureLayer_management(in_feat4,out_layer4)
                print "Make Feature Layer for " + fc

                #make feat layer of fc_endswith_d3
                in_feat3 = buffs +"\\"+ fc[:-2]+"d3"
                out_layer3 = fc[:-2]+"d3"
                arcpy.MakeFeatureLayer_management(in_feat3, out_layer3)
                print "Make Feature Layer for " + in_feat3
### ERROR IN THIS SECTION ####
                for ID4 in APR16_ID4:
                    where_clause4=""" "APR16_ID" = """ + "'"+ID4+"'"
                    arcpy.SelectLayerByAttribute_management(out_layer4,"NEW_SELECTION",where_clause4)
                    arcpy.SelectLayerByAttribute_management(out_layer3,"NEW_SELECTION",where_clause4)
                    erase_out= "C:\\BOEM\\Donuts\\Donuts_d4_output.gdb\\"+"don_d4_"+ID4
                    ### Error in the Erase function here:
                    arcpy.Erase_analysis(out_layer4, out_layer3,erase_out,"")
                    print "Erase "+ out_layer3 +" from "+out_layer4+" to produce "+erase_out
                    print arcpy.GetMessages()
                    print "___"


texts: Importing QGIS layer to SpatiaLite table

texts:
I have a Spatialite database with a table called spatialite_objects, with spatial and non-spatial data.

Another layer is an Oracle Spatial Layer (id, geom) with 2 table joins (based on 2 csv files) called oracle_objects.

I use the DB Manager to import this layer to my Spatialite database, successfully.

Now, I want to do this with a QGIS Plugin using Python.

I know to create tables and doing queries with the pyspatialite library (from pyspatialite import dbapi2 as db). I think about to query the Oracle layer and to fill up the Spatialite table.

Is this the same way the DB Manager does the import or are there more "direct" commands to import a layer from the Canvas directly to a Spatialite database?


texts:
I want to add custom projection to the GDAL library version 3.
In earlier versions (with Ubuntu 18.04) we directly added the definition to the /usr/share/proj/epsg file
In GDAL 3 the files do not exist any more, how can we add the custom projection?


texts: From QGIS 3.4.5 and 3.6 asks password repeadly to connect to PostgreSQL 10

texts:
Action: At first time when opening project on QGIS 3.4.5 and 3.6 versions needs credentials to connect to PostgreSQL DB.

Problem: But after that when connection is success with correct credentials the qgis repeatedly asks the same credentials again to connect to the same DB when the maps is zooming and moving and etc.

With QGIS 2.18, 3.4.0, 3.4.1 versions all works fine and don't asking repeatedly credentials after first correct time.

Maybe someone knows whats is wrong?


texts: QGIS creating a fixed distance buffer over selected areas not the entire layer

texts:
what im doing is making a fixed distance buffer and it's working splendidly but i need the buffer to be created over certain selected areas not the entire layer, and at the creation wizard is only offers the the option of selecting the entire layer, and i need some areas within that layer to have a different buffer and some none at all.

is there a way around that?


texts: How to set the extent with gdal_rasterize, the "-te" does not work

texts:
I am using gdal_rasterize to clip mask by shapefile. However, the parameter "-te" does not work.

I use subprocess to call GDAL util:

import subprocess

# define paths to raster and vector
inraster = r'I:\HY_project\HY1B-JPG2\2011\deepwater\test\H1B_OPER_OCT_L1A_20110120T_011319707_geo.img'
inshape = r'I:\HY_project\HY1B-JPG2\2011\deepwater\test\H1B_OPER_OCT_L1A_20110120T_011319707_geo.shp'

# the extent of shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shp, 0)
layer = dataSource.GetLayer()
feature = layer[0]
geom = feature.GetGeometryRef()
extent = geom.GetEnvelope()
xmin = extent[0]
xmax = extent[1]
ymin = extent[2]
ymax = extent[3]

cmd = 'gdal_rasterize -burn -999 -te '+str(xmin)[0:7]+' '+str(ymin)[0:7]+' '+str(xmax)[0:7]+' '+\
     str(ymax)[0:7]+' '+inshape + ' ' + inraster

cmd = 'gdal_rasterize -burn -999 -i -te 130.827 24.8857 140.357 26.8057 '+ inshape + ' ' + inraster

# run command
subprocess.call(cmd, shell=True)

If the parameter "-te" is not set then the following command line works:

cmd = 'gdal_rasterize -burn -999 -i '+ inshape + ' ' + inraster


texts: Clipping raster with vector layer using MapInfo?

texts:
With MapInfo 12, i would like cut my raster layer .tab with one of my vector layer.

How to do this basic task with MapInfo?

Is it possible with mapinfo to convert .tab raster table toward any standard formats (ecw, tif, etc)?


texts:
I want to use the result of OD cost matrix for my further calculations. However, whne I try to extract the sublayer "lines" it returnes an error

AttributeError: 'Layer' object has no attribute 'listLayers'

However a line before the method saveACopy (file_name) worked on the same layer file.

import arcpy
import os
def execute(self, parameters, messages):
    # Check out the Network Analyst extension license
    arcpy.CheckOutExtension("Network")

    # Set environment
    output_dir = r"D:\GISworkspace"
    arcpy.env.workspace = os.path.join(output_dir,'DataTest.gdb')  # Setting the database, where the output table will be created
    out_path = arcpy.env.workspace                                 # Variable so that it is easier to operate with the path

    #***************************************************************************

    #Set local variables
    layer_name = "ODcostMatrix"
    impedance = "Length"
    output_layer_file = os.path.join(output_dir, layer_name + ".lyrx")
    network = parameters[0].value
    nodes   = parameters[1].value

    result_object = arcpy.na.MakeODCostMatrixLayer(network, layer_name, impedance)

    #Get the layer object from the result object. The OD cost matrix layer can
    #now be referenced using the layer object.
    layer_object = result_object.getOutput(0)

    #Get the names of all the sublayers within the OD cost matrix layer.
    sublayer_names = arcpy.na.GetNAClassNames(layer_object)
    #Stores the layer names that we will use later
    origins_layer_name = sublayer_names["Origins"]
    destinations_layer_name = sublayer_names["Destinations"]

    #Load the BS locations (nodes) as both origins and destinations.
    arcpy.na.AddLocations(layer_object, origins_layer_name, nodes)
    arcpy.na.AddLocations(layer_object, destinations_layer_name,nodes)

    #Solve the OD cost matrix layer
    arcpy.na.Solve(layer_object)

    #Save the solved OD cost matrix layer as a layer file on disk
    layer_object.saveACopy(output_layer_file)

    #Get the Lines Sublayer (all the distances) TODO fix the Atrribute Error: no method listLayers
layer_object.listLayers(sublayer_names["ODLines"])[0]


texts: Isolate the largest patch in a raster in GRASS?

texts: I have a raster with a series of patches (contiguous areas with same value). How do I isolate the largest patch? I know I could run r.stats, look through the long list of calculated areas, figure out which one is the largest and then use r.reclass to isolate it, but I'm sure there has to be an easier or more efficient way to do this.


texts: How to digitize rectangles with snapping in QGIS?

texts:
Trying to digitize a floor plan for a building (just need polygons for each room) in QGIS (shapefile layer).  I started just trying to draw the features, but want to align polygons so that there is no gap or overlap.   That led me to the Snapping features which helped a lot.  But having trouble making sure my rectangles have 90 degree corners.

Then I discovered the rectangles ovals digitizing tool which solves that problem.  Unfortunately--it seems to ignore the "snapping" that I got when just drawing features node to node.  Is there a way to get both snapping and rectangles in QGIS?


texts: Why does Modify Vertices seem to fail for no reason in ArcGIS Pro?

texts:
I could be working on a project for a while, between two maps and having no trouble modifying features. Then suddenly I will try to adjust the vertices on a line, and I get an error message:

  modify vertices failed

Then, I can't seem to modify any features! I've tried restarting the program, and my computer. Has anyone else encountered this before? I have tried scouring the forums but couldn't find any other reports of this error. I'm working in ArcGIS Pro 2.3.0 on Windows 10. Here's a screen shot of my error and the image of the feature in my project that I am trying to modify.


texts: How to drape image on surface using open source software?

texts:
I have a image with a known projection - say an aerial photo. I also have a raster surface grid with the same projection - say a DEM. Both are in the same projection. Both have the same extent.

How do I "drape" the image on top of the surface using open source software?

The end result would be an oblique-like image that you typically see when you use google earth (or similar 3D earth-observation software).

This question does what I want to do (using Esri), but get better results - What can be used to drape an High res Image over a DEM and get high quality results? And this is a similar (but different) question for QGIS -How to drape a simple jpeg image over dem of some area?

There doesn't seem to be a way to do it in QGIS. I think GRASS should be able to do it, but I haven't figured out a way by looking around or googling. ParaView seems like a possibility, but overly complicated.


texts: Running GRASS GIS on Text-based interface alone

texts:
I'm required to run GRASS GIS functions on a CentOS 7 server using the GRASS Text-based interface as it is not possible to use the GUI. I basically want to import some shapefiles as vector layers and run a few network analysis functions. I wanted to see if the v.import command works but I just ended up with
bash: v.import: command not found

I'm fairly new to this so I'm not quite sure what I need to do. Is there no way to perform functions on a layer using the text based UI (CLI)


texts: How can I split one band from 3 stacked layer image by ERDAS IMAGINE?

texts:
I have a 3 layers stacked image, each layer represents a band ratio for lithology extractions for geological purposes. I need to separate each layer separately
how can I do this by using erdas imagine


texts: Which rule to follow for defining an intersection as an ordered pair of streets?

texts:
Explanation

Sometimes, for the sake of simplicity or speed, I need to list a location as a pair of intersecting streets. Let's use Luckie Street Grocery Store in Atlanta, Georgia as the example for now. I might list it as being at Luckie St NW & Mills St NW.

Often, I wonder which street should be listed first and why. There are some heuristics I could use to answer this question, such as:

Cartesian coodinates — Luckie St NW & Mills St NW

Using this, I would list the vertical street first, followed by the horizontal street.

Geographic coordinates — Mills St NW & Luckie St NW

Using this, I would list the horizontal street first, followed by the vertical street. (There is considerable debate in the scientific community about this standard, but that is off-topic.)

"Better known" steet first — Luckie St NW & Mills St NW

It would depend on how well one understands the locale to successfully utilize this method. I don't live in Atlanta or know that area well, but I took a stab at it anyway using zoomed out Google Maps as my reference.

Other ideas

I could use alphabetical order (by street name or by type of street), but I could encounter conflict by following this rule (two streets with the same name, two streets of the same type).

Additionally, there are plenty of intersections which are not closely aligned to North-South (longitude) and East-West (latitude) lines (i.e. they are more in the shape of an X on a map). An example of this in the same city is Poplar St NW & Broad St NW. This presents a problem for the first two methods I mentioned.

Question

Including a link to one or more authoritative sources, what is the preferred method of determining the order of streets when listing an intersection (and, most importantly, why)?


texts: How to calc area using OpenLayers with UTM coordinates


texts: Spliting the column in CSV its showing attribute Error in QGIS plugin

texts:
I am trying create a new plugin in QGIS. I have my own script its work in normal Jupyter Notebook, then I planned to convert into as QGIS Plugin.The dout/error occur while spliting the columns.
    def run(self):
    """Run method that performs all the real work"""
    if self.first_start == True:
        self.first_start = False
        self.dlg = DegreesminutessecondstodecimaldegreesDialog()
        self.dlg.pushButton.clicked.connect(self.select_output_file)

    # Fetch the currently loaded layers
    layers = QgsProject.instance().layerTreeRoot().children()
    self.dlg.MapLayerComboBox.addItems([layer.name() for layer in layers])
#Connecting to selected layer
    def layer_changed(layer):
        self.dlg.FieldComboBox.setLayer(layer)
        self.dlg.FieldComboBox_2.setLayer(layer)
    self.dlg.MapLayerComboBox.layerChanged.connect(layer_changed)

    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed

    if result:
        def D_D():
            DD = self.dlg.MapLayerComboBox.layerChanged.connect(layer_changed)

            lat = self.dlg.FieldComboBox.currentField()
            lon = self.dlg.FieldComboBox_2.currentField()

            g = pd.DataFrame(lat.str.split(expand=True).tolist(),
                             columns=['1', 'D1', 'D2', '4', 'M1', 'M2', '7', 'S1', 'S2', '10'])

I run its show error in
g = pd.DataFrame(lat.str.split(expand=True).tolist(),
AttributeError: 'str' object has no attribute 'str'

When I remove str before split
g = pd.DataFrame(lat.split("").tolist(),
ValueError: empty separator

Can you give me some ideas to make this work?


texts: Enhance Vegetation Index ( EVI ) equation in ENVI IDL

texts: I'm using ENVI 5.3 to find EVI, when I'm applying the equation it doesn't run, anybody can provide me the exact equation or formula for running the index. For more clarification, I'm using Landsat 8 OLI.


texts: RAT and cell values in R's raster-objects

texts:
I really struggle to understand some basic concepts in the structure of how R's raster-objects save values. Namely the RasterLayer, the RasterBrick and the RasterStack.

I understand that the values (when manually creating a raster) can just be of type numeric, integer, logical or factor. When I create a raster with numeric values like this:

raster_numeric = raster(nrows = 6, ncols = 6, res = 0.5,
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = seq(0.1,3.6,0.1))

I get this structure:

class      : RasterLayer
dimensions : 6, 6, 36  (nrow, ncol, ncell)
resolution : 0.5, 0.5  (x, y)
extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 0.1, 3.6  (min, max)

And apparently no attributes-slot.

Whereas, when I use factorial values (copied from the amazing Lovelace et al. book) I get, what to my understanding is a Raster-Attribute-Table (RAT). So why is this?

class      : RasterLayer
dimensions : 6, 6, 36  (nrow, ncol, ncell)
resolution : 0.5, 0.5  (x, y)
extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 1, 3  (min, max)
attributes :
 ID VALUE
  1  clay
  2  silt
  3  sand

What exactly is this RAT and can raster-objects really save just one value because in the end each cell can be references to a row in the RAT and thereby hold many values of different type?


texts: Why does the map take so long to display and appear distorted?

texts:
I have map rendered from geoserver but it takes too much time, again shown map is half-shown and distorted as the following image, SOMETIMES THE MAP DOES NOT EVEN DISPLAY

 My JAI setting is as below

<jai>
    <allowInterpolation>false</allowInterpolation>
    <recycling>false</recycling>
    <tilePriority>5</tilePriority>
    <tileThreads>7</tileThreads>
    <memoryCapacity>0.5</memoryCapacity>
    <memoryThreshold>0.75</memoryThreshold>
    <imageIOCache>1</imageIOCache>
    <pngAcceleration>true</pngAcceleration>
    <jpegAcceleration>true</jpegAcceleration>
    <allowNativeMosaic>false</allowNativeMosaic>
</jai>

URL IS

http://localhost:8080/geoserver/it.geosolutions/wms?service=WMS&version=1.1.0&request=GetMap&layers=it.geosolutions:thanaboundary&styles=&bbox=88.1099243164063,21.1824684143066,92.5484466552734,26.4899482727051&width=428&height=512&srs=EPSG:4326&format=application/openlayers

My PostGIS DataBase Content in SQL is as below

  -- Table: thana

  -- DROP TABLE thana;

  CREATE TABLE thana
  (
        gid serial NOT NULL,
        objectid integer,
        recid numeric,
        geocode numeric,
        div_code numeric,
        div_name character varying(254),
        dist_code numeric,
        dist_name character varying(254),
        thana_name character varying(254),
        thana_code numeric,
        shape_leng numeric,
        shape_area numeric,
        noofbed double precision,
        bedoccupen double precision,
        dailyavgop double precision,
        geom geometry(MultiPolygon,4326),
        CONSTRAINT thana_pkey PRIMARY KEY (gid)
  )
 WITH (
       OIDS=FALSE
      );
       ALTER TABLE thana
       OWNER TO postgres;

  -- Index: thana_geom_gist

  -- DROP INDEX thana_geom_gist;

          CREATE INDEX thana_geom_gist
          ON thana
          USING gist
          (geom);

My geoserver web page

Would any one with any tweak please


texts: How to visualize choropleth with temporal component in browser

texts:
I have a vector dataset with monthly water risk data per river basin. The data is available for each month (January - December) and I would like to visualize this in a browser. The geometries are constant (river basins do not change).

The  data is  structured as two files:  a shapefile with unique identifiers and a .csv file with the unique identifiers, month and value.

The preferred way to display one month of data is using choropleths. I think a very powerful way to visualize the temporal component is to have a time slide and "play" button.

I've looked at Torque JS, Mapbox GL JS however the examples do not include choropleth examples. I can imagine other people might've solved this challenge already.


texts: Mapbox Vector Tiles appear to lack accurate intersection nodes?

texts:
Mapbox Vector Tiles appear to lack accurate intersection nodes.

Is this a bug?

i.e. should the tile's geojson include a shared node between two intersecting streets?

Details:
Using MVT's from http://a.tiles.mapbox.com/v4/mapbox.mapbox-streets-v8 and converting them to json via MB's https://github.com/mapbox/vt2geojson, the street data lacks many intersections that are correctly included in the original OSM Overpass data.

By this I mean that there should be a node at the intersection that is shared by both street's LineString's.

Is this a bug, either in the spec or in the MB implementations?


texts: Which basename i am going to use as prefix in i.landsat8.swlst in GRASS GIS?

texts:
I want to use i.landsat8.swlst add-on in GRASS GIS to calculate LST from Landsat8 data for my study area.

In the following command:

i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -n

Which Basename i am going to provide as prefix for landsat8 data?


texts: R which census division a longitude, latitude pair falls into

texts:
I have the dataframe below

> head(cities)
  location_id longitude latitude
2       1  -79.3942  43.6724
3       2 -123.1070  49.2540
4       3  -71.0589  42.3601
5       4  -84.3880  33.7490
6       5  -97.7431  30.2672
7       6  -76.6122  39.2904

I downloaded the census division boundary shapefile from statcan which I read using

shapefile <- readOGR(".", "gcd_000a11a_e")

http://www.statcan.gc.ca/access_acces/alternative_alternatif.action?l=eng&dispext=zip&teng=gcd_000a11a_e.zip&k=%20%20%20%2014372&loc=http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/gcd_000a11a_e.zip

I'm trying to find the census division each pair of coordinates fall into.

> summary(shapefile)
Object of class SpatialPolygonsDataFrame
Coordinates:
         min       max
x -141.01807 -52.58230
y   41.67695  89.99943
Is projected: FALSE
proj4string :
[+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0]
Data attributes:
    CDUID              CDNAME             CDTYPE             PRUID              PRNAME
 Length:293         Length:293         Length:293         Length:293         Length:293
 Class :character   Class :character   Class :character   Class :character   Class :character
 Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character
>


texts: Export Map Error from OpenLayers 3 in Chrome?

texts:
Im using this js file to export map in Openlayers 3 in Google Chrome but i'm getting en error saying that Uncaught DOMException: Failed to execute 'toDataURL' on 'HTMLCanvasElement': Tainted canvases may not be exported. I also set an cross Origin attribute as anonymous. The function of export button is as follow

var exportPNGElement = document.getElementById('export-png');
exportPNGElement.crossOrigin = "anonymous";
    if ('download' in exportPNGElement) {
      exportPNGElement.addEventListener('click', function(e) {
        map.once('postcompose', function(event) {
          var canvas = event.context.canvas;
          exportPNGElement.href = canvas.toDataURL('image/png');
        });
        map.renderSync();
      }, false);
    } else {
      var info = document.getElementById('no-download');
      /**
       * display error message
       */
      info.style.display = '';
    }


texts:
I am catching the selectionChanged signal and passing the selected layer to the registered method. Once here, I want to know the features that were selected. I do this with a call to layer.selectedFeatures().

The return type from this call is a QgsFeatureList. As a test, I'm trying to print the values of the features in the QgsFeatureList that was selected, or just display them in a simple way. I find that neither QgsFeatureList nor QgsFeature has any functions returning a QString.

How can I print the values of the features that were selected?


texts:
I have developed a fishnet grid for my study area.

I have various datasets for the area and want to allocate a value (from one of the pre-determined datasets) to each polygon within the Fishnet.

How do i display that the data from one source falls into the various polygons?


texts: leaflet tilelayer from localhost geoserver

texts:
I am attempting to construct a tilelayer link from a local geoserver. the layer is called county and I seem to have set it up correctly to be called as a tilelayer

now as I am following the leaflet documentation http://leafletjs.com/reference-1.2.0.html#tilelayer

I am getting confused on how to construct the URL

here is my code

  <title>A Leaflet map!</title>
  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>

nothing shows up when I display this in google chrome.

in the console I get this error

how do I properly construct this tilelayer URL?


texts: Pirmary key constraint on objectid field

texts:
I am trying to copy data from a PostgreSQL (9.4) DB table to an SDE table. Both have an objectid field. If I copy the data from the source table with the OID field, I get an error because the values in OID in the destination table have duplicates. If I copy the data without the OID field, I get an error for not respecting the not null constraint on the destination OID.

I tried to remove the primary key constraint and the not null constraint on the OID with an alter table and even if through PGadmin I can see that there is no pkey constraint, as I run my script, I still get the error saying that the not null constraint can't be removed because objectid is in a primary key.

I would like to avoid reassigning a new OID to the entries with a cursor in the source table, as there is a lot of them. Is there any way I could bypass those constraints?


texts: Seeking syntax for arcpy.mapping.Layer() in ArcPy for ArcGIS Pro?

texts:
I'm working on a project where I am porting Python Tools created for ArcMap Desktop into ArcMap Pro but I am encountering an issue when I get to lines of the original code that read:

insertLayer = arcpy.mapping.Layer(input), which aims to create a Layer object to manipulate/process.

Having presumed product compatibility, I tried insertLayer = arcpy.mp.Layer(input) but that didn't work - I get the error message module 'arcpy.mp' has no attribute 'Layer'

From the help pages, I can see there are two ways to access Layer object:

Using listLayers method of Map object, which is dependent on an
existing project file (.aprx) with the map having the target layer
Using MakeFeatureLayer_management, which creates a Layer from a feature class and let you to use of a string (i.e., the second argument) in place of a layer object but its use is generally limited to Select from Layer operations.

What is the best way to port the Layer object, which is created in ArcGIS Desktop, to Pro?


texts: trim a global quadtree by population density?

texts:
I'm looking for a way to match people so as to find other members that are 'nearby' - not for those dreadful 'girls near you' ads, don't worry ;_)

Of course 'nearby' in the country is a different thing than when in the city.

I'm thinking that an efficient way to do this would be to use a quadtree for their location, but use deeper level divisions when the population density is higher. In other words try to 'balance' the tree so that every leaf node has roughly the same population within it.

I was thinking that someone may have already done this but can't find such a thing on the web.

If not, is there an alternative approach I should look into?


texts: JavaScript to Python, how to translate bitewiseAnd?

texts:
I'm using the Earth Engine API in a Jupyter notebook and I try to translate from JavaScript to Python a cloud masking function (the one given with the dataset)
/**
 * Function to mask clouds based on the pixel_qa band of Landsat SR data.
 * @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};

I already fell in the trap of and, or and not that are protected in Python. I'm now facing an error with the bitewiseAnd keyword:
def clouMask_landsat(image):
    """
    Function to mask clouds based on the pixel_qa band of Landsat SR data.

    Args:
        image(ee.Image): image Input Landsat SR image
    Returns:
        (ee.Image): Cloudmasked Landsat image
    """
    qa = image.select('pixel_qa')

    # If the cloud bit (5) is set and the cloud confidence (7) is high
    # or the cloud shadow bit is set (3), then it's a bad pixel.
    cloud = qa.bitewiseAnd(1 << 5).And(qa.bitewiseAnd(1 << 7)).Or(qa.bitewiseAnd(1 << 3))

    # Remove edge pixels that don't occur in all bands
    mask2 = image.mask().reduce(ee.Reducer.min())

    return image.updateMask(cloud.Not()).updateMask(mask2)

AttributeError: 'Image' object has no attribute 'bitewiseAnd'

Is it simply a & or is it something else entirely?


texts: WFS GetFeature request with box failed on OGR datastore

texts:
I prepared an OGR datastore and imported features from a GeoJSON data set using REST.

If I now send a wfs?request=GetFeature&bbox=....., I get an error message:

java.lang.NullPointerException: PropertyDescriptor is null
did you request a property that does not exist?

I traced it down to: ComplexTypeImpl:70 where an invalid property occurs.

This happens since an additional bogus attribute with name="" was added at: OGRFeatureSource:213 Set<String> extraAttributeSet = ...

Those extra attribute is extracted from the query filter which is a BBOXImpl. I guess the intention is to get "the_geom" for the filter envelope to apply on.

The BBOX of the query is generated by a geoserver class at:
BaseFeatureKvpRequestReader:354

Here a comment suggests:

// JD: use "" so that it applies to all geometries
String name = "";

Seems that other data stores have no problem with this.

Any idea?


texts: How to make maps move when you create a Google map using OpenLayer 2

texts:
<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=0">
    <meta name="apple-mobile-web-app-capable" content="yes">
    <title>OpenLayers Google Layer Example</title>
    <link rel="stylesheet" href="../theme/default/style.css" type="text/css">
    <link rel="stylesheet" href="../theme/default/google.css" type="text/css">
    <link rel="stylesheet" href="style.css" type="text/css">

Google Layer Example

        Google

        Demonstrate use of the various types of Google layers.

        For best performance, you must be using a version of the Google Maps
        API which is v2.93 or higher. In order to use this version of the API,
        it is best to simply set your application to use the string "v=2" in
        the request, rather than tying your application to an explicit version.

       In order to position the Google attribution div in the default location,
       you must include the extra theme/default/google.css stylesheet.

You can zoom in and out. But moving is impossible. I don't know how to fix the code.

Written with Google Translator, it may be strange.


texts:
I have a dataset of road lines obtained from digitized aerial imagery and I want to automate the generalization procedures using arcGIS model builder. I am using the ArcGIS 10.3 desktop suite. I have come across the generalization toolset in ArcGIS. But I'm not sure as I may not have understood the theory in generalizing roads.

What procedures can I use in a workflow?

And especially, since I'm generalizing from my source data at 1 in 2500, to an output of 1 in 20k how can i calculate the parameters to use for a 20k map generalization module?


texts:
I downloaded an earthquake shake map for USGS with KMZ file.

It is appears in Google Earth as below.

I would like to use this intensity in ArcMap. How can I do this?


texts: Making labels stay at some defined position in QGIS?

texts:
I am using QGIS 3.2 Bonn.

If I add a label (and I have a lot of labels, for every field, river, path, etc...) and if I zoom out all labels are crossing each other and it is hard to know what is happening then. In the first picture it is like a big mess caused when I zoom out, and the 2nd picture is an example of how I want it to stay forever.

Is there a way that labels will not "stick out" when I zoom out of the map, and will just stay until I zoom back in?


texts: Convert OL3 drawn polygon to WKT for inserting to PostGIS

texts:
Essentially what I am trying to accomplish is to have a map where the user can draw a geometry, fill in a form and have the record, along with the geometry, saved to PostGIS.  I'm new to Javascript (and newer still with openlayers3) so everything I've managed to piece together is from random corners of the internet.

Below is my full script, which currently can get a map, have a geometry drawn on it and finally show the coordinates in an empty < div id="wktgeom" > from a listener on 'drawend'. I've tried this solution:

var ol3Geom = vectorSeleccion[x].getGeometry();
var format = new ol.format.WKT();
var wktRepresenation  = format.writeGeometry(ol3Geom);

without success.  I've tried other things but i get an output to the < div > as [object Object].  This could easily be my lack of knowledge and I'm putting code in the wrong place (currently i'm sticking everything into the listener which I can get the coordinates to output from).

Can someone please show me how to fix my code so that a single polygon drawn on the map outputs the geometries coordinates as WKT (it will only ever be a single geometry), and also any answers have some form of context to my script so I 100% know I'm putting it in the right place?

    <script type="text/javascript">
      var source = new ol.source.Vector({wrapX: false});

      var vector = new ol.layer.Vector({
        source: source,
        style: new ol.style.Style({
          fill: new ol.style.Fill({
            color: 'rgba(255, 155, 155, 0.4)'
          }),
          stroke: new ol.style.Stroke({
            color: '#ff0000',
            width: 2
          })
        })
      });

      var map = new ol.Map({
        target: 'map',
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          }), vector
        ],
        view: new ol.View({
          center: ol.proj.fromLonLat([-0.7594, 52.0406]),
          zoom: 8
        })
      });

      var typeSelect = document.getElementById('type');

  var draw; // global so we can remove it later
  function addInteraction() {
    var value = typeSelect.value;
    if (value !== 'None') {
      var geometryFunction, maxPoints;
      if (value === 'Square') {
        value = 'Circle';
        geometryFunction = ol.interaction.Draw.createRegularPolygon(4);
      } else if (value === 'Box') {
        value = 'LineString';
        maxPoints = 2;
        geometryFunction = function(coordinates, geometry) {
          if (!geometry) {
            geometry = new ol.geom.Polygon(null);
          }
          var start = coordinates[0];
          var end = coordinates[1];
          geometry.setCoordinates([
            [start, [start[0], end[1]], end, [end[0], start[1]], start]
          ]);
          return geometry;
        };
      }
      draw = new ol.interaction.Draw({
        source: source,
        type: /** @type {ol.geom.GeometryType} */ (value),
        geometryFunction: geometryFunction,
        maxPoints: maxPoints
      });

      //Here is where I've tried outputting WKT using the other solution without success
      draw.on('drawend',function(e){
        document.getElementById('wktgeom').innerHTML = e.feature.getGeometry().getCoordinates();
      })
      map.addInteraction(draw);
    }

  }

  /**
   * Handle change event.
   */
  typeSelect.onchange = function() {
    map.removeInteraction(draw);
    addInteraction();
  };

  addInteraction();
</script>


texts: Converting none to zero (0) in geodatabase table using ArcPy?

texts:
I have an attribute table looks like this

How do I convert Nulls into 0?

I've tried to set this table non nullable type but got this:

while the official doc says this:

At least one of them is wrong.

Then I tried many ways to do the job, like this:

RValues=arcpy.SearchCursor(sr)

for RValue in RValues:
    #print (RValue.getValue(field))
    if RValue.getValue(field) == None:
       #print (RValue.getValue(field))
       return 0
    else
       return
       float(RValue.getValue(field))

print(RValue.getValue(field))

This gives error 999999, don't know why, possible a bug.

Then I tried this:

RValues=arcpy.SearchCursor(sr)
print (RValues)

for RValue in RValues:
    if RValue.getValue(field) == None:
        print (RValue.getValue(field))
        arcpy.CalculateField_management(sr, "RValue",  esp, "PYTHON3")
    print (RValue.getValue(field))

Well, the first print gives a

<geoprocessing cursor object object at 0x0000025F5C161E70>

the second print does isolate all none items, as there is no null in python.

Here is the trick:

the third print either gives error or turn everything into null or 0,focus on everything, expressions I use are as follows:

esp="!RValue!.replace(None, 0)"

esp2="!RValue! * 0"

esp3="set field = 0 where field is null"

esp4="NVL(RValue,0)"

esp5="select nvl(RValue, 0) from sr.RValue"

esp6="None"

How is that even possible when the if before it already isolate only those nones?