texts:
I am trying to ultimately compare the gate height to the corresponding fence height located nearby. So far I have come up with the following:

gates = "gate_line"
fence = "fence_line"

arcpy.Near_analysis(gates, fence)

with arcpy.da.SearchCursor(gates, ["NEAR_FID", "gate_ht"]) as cursor1:
  with arcpy.da.SearchCursor(fence, ["OBJECTID", "fence_ht"]) as cursor2:
    while row in cursor1:
        for rows in cursor2:
             print "Fence Height: {0}, Gate Height: {1}".format(rows[1], row[1])

Obviously, I am not to a comparison stage yet, still trying to get the output to read properly to know that I am seeing results from both the feature classes represented. With this code block the fence heights are displaying correctly, however the gate heights are all the same height (they are not in the data).

Suggestions?

Edit:

This is what I have ended up with and it is working:

with arcpy.da.UpdateCursor(gates, ["NEAR_FID", "gate_ht", "OBJECTID", "SAME_HEIGHT"]) as cursor1:
    for row in cursor1:
        sqlStr = '"OBJECTID" =' + str(cursor1[0])
        with arcpy.da.SearchCursor(fence, ["OBJECTID", "fence_ht"], sqlStr)  as cursor2:
            for rows in cursor2:
                if row[1] != rows[1]:
                    row[3] = "False"
                cursor1.updateRow(row)


texts: Is there any place in which I can find a list of invalid database characters for ArcGIS 10?

texts:
Relating to my question about the split tool here, as suggested I ask a new question: is there any place with a list of invalid characters to use in a database in ArcGIS 10?

edit: It is in relation to the error I got before which is "The name starts with an invalid character". It didn't say where beyond that.


texts: Add more info in the attribute table after executing "Line intersection" in QGIS

texts:
I have executed the "Line intersection" tool on a road shapefile. This is under Vector > Analysis Tools > Line intersections. I get an output with attribute table having just two columns: type_1 and type_2.

Is it possible to get more info/attributes out of this? For example, the node_ids of the intersections and perhaps their location (latitude and longitude)?


texts: Force an ArcGIS API for JavaScript WMSLayer to use EPSG: 900913 (instead of 102100)

texts:
I have a use case where an external WMS server can serve maps in the Web Mercator projection using EPSG code 900913, but won't recognize ESRI's 102100 code.

My question is: can I force all WMS requests to use EPSG code 900913?

To illustrate, I have a simple map page:

<!DOCTYPE html>

        <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
        <title>Test</title>

        <link rel="stylesheet" href="http://js.arcgis.com/3.8/js/esri/css/esri.css">

Where test.js creates a simple map and adds a WMSLayer to it:

var wmsAddress = "http://some.server.com/wms";
var map, srs, extent;

require(["esri/map",
         "esri/layers/WMSLayer",
         "esri/layers/WMSLayerInfo",
         "esri/geometry/Extent",
         "esri/geometry/Point",
         "esri/SpatialReference",
         "dojo/domReady!"],
        function (Map, WMSLayer, WMSLayerInfo, Extent, Point, SpatialReference) {

    srs = new SpatialReference({wkid:900913});
    extent = new Extent(
            1362550.0288221193, 7616458.979152972,
            2281987.5791366613, 9525359.127428867,
            srs);
    map = new Map("mapDiv", {
        center: new Point(1800000, 7900000, srs),
        extent: extent
    });

    var wmsService = new WMSLayer(wmsAddress, {
        resourceInfo: {
            extent: extent,
            layerInfos: [new WMSLayerInfo({name:"layerName",title:"TestMap"})]
        },
        visibleLayers: ["layerName"]
    });
    map.addLayers([wmsService]);
});

When loading the application, the request that is sent to the WMS server looks like this (viewed in a browser console):

  http//some.server.com/wms?SERVICE=WMS&REQUEST=GetMap&FORMAT=image/png&TRANSPARENT=TRUE&STYLES=&VERSION=1.3.0&LAYERS=layerName&WIDTH=1904&HEIGHT=400&CRS=EPSG:102100&BBOX=-2743182.35289663,6945549.925862053,6343182.35289663,8854450.074137948

Which results in this response from the server:

<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://some.server.com:80/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">
    <ServiceException code="InvalidCRS">
        Error occurred decoding the espg code urn:x-ogc:def:crs:EPSG:102100 No code "EPSG:102100" from authority "European Petroleum Survey Group" found for object of type "IdentifiedObject".
    </ServiceException>
</ServiceExceptionReport>


texts: Converting a rasterfile to a point shapefile with all columns from the rasterfile

texts:
I'm trying to convert a raster with multiple columns to a point shapefile to get rid of the "count" column but unfortunately the tool "raster to point" works just for one column so I tried to modify this code to write the atribute table to a csv file but it still doesn't work.

I just get a lot of zeros in my outfile. Why am I get just a bunch of zeros?

I'm using ArcMap.

This is the modified code:

outfile = Workspace + "/table"

fOut = open(outfile, "wb")
fOut.write('ID' + ',' + 'Flowaccumaltion' + ',' + 'Flowlength\n')
rstArray = arcpy.RasterToNumPyArray(Combined_Raster)
rows, cols = rstArray.shape
for rowNum in xrange(rows):
    for colNum in xrange(cols):
        if colNum == 1 or  colNum == 3:
            value = rstArray.item(rowNum, colNum)
            fOut.write(str(value) + ',')
        elif colNum ==4:
            value = rstArray.item(rowNum, colNum)
            fOut.write(str(value) + '\n')
fOut.close()

This is the raster attribute table I want to convert

And I would like to get something like this:

ID  , Flowaccumulation,  Flowlength
57  ,    727          ,      1
58  ,     727         ,      1
 0  ,   76073         ,   1008
 1  ,   14919         ,    622
 .  ,      .          ,      .

But if there is a tool which can convert a raster to a point shapefile with all fields except for the "count" field then this would be even better.


texts: Understanding ArcGIS Server connection types

texts:
Being relatively new to ESRI ArcGIS Server I have a question regarding the storage of ESRI SDE Databases - AGS 10.0 Enterprise and ArcSDE 9.3.

Firstly the naming of the whole structure gets confusing. SDE as far as I know, is only the software which "spatially enables" a Database (in my case MS SQL 2008 R2). So when I create a geodatabase in MSSQL I have to do this through the post install program. Only then will I have the ability to view/update/delete in the geodatabase either with ArcEditor/Info.

Im my "inherited" system that I have to now continue using, I have some data coming from MSSQL Server and some from Geodatabases which are on another PC.

Question 1: What type of connection is that on the pc? Is that a so called "direct" sde connection?

Question 2: Some other colleagues have told me, that Rasters are delivered quicker in AGS Services if they are not in the MSSQL DB. Is this true?

Thanks for any help,

Robert


texts: How to iterate nearest polygons until threshold is reached?

texts:
I have a point layer and a polygon layer (representing buildings). I created a buffer around the points and extracted those polygons which intersect the buffer.

The points contain an attribute with a certain value (e.g. 100). The polygons also contain an attribute with varying values (e.g. 5.5, 7.9 etc.).

What I would like to do is to iterate the intersected polygons in terms of distance to the point (i.e. the polygon closest to the point would be first, the polygon second closest would be second etc.). And each time, the polygon value would subtract the point value until it can no longer be subtracted (it cannot be less than 0) in which case it should stop iterating.

If possible, it would also be ideal where if a polygon has been iterated over, a new field would contain a boolean value 1 and all those which have not been iterated over contain 0.

I'm hoping the result would look similar to this:

How could this be achieved by either SQL (virtual layer) or PyQGIS?

The following are pieces of code which only extracts the intersected polygons at 200m from the points:

SQL

SELECT * FROM
Buildings d,
point c
WHERE ST_INTERSECTS(d.geometry, ST_BUFFER(c.geometry,200))

PyQGIS

def func(point_layer, polygon_layer, spacing):
    # Create buffer
    buffer_layer = QgsVectorLayer("Polygon?crs=epsg:27700", 'Buffers' , "memory")
    with edit(buffer_layer):
        buffer_layer.addAttribute(QgsField('ID', QVariant.Int))
        feat_list = []
        for f in point_layer.getFeatures():
            poly = QgsFeature()
            f_buffer = f.geometry().buffer((spacing), 99)
            poly.setAttributes([1])
            poly.setGeometry(f_buffer)
            feat_list.append(poly)
        buffer_layer.addFeatures(feat_list)
    buffer_layer.updateExtents()
    buffer_layer.updateFields()
    QgsMapLayerRegistry.instance().addMapLayers(buffer_layer)
    # Intersect with polygon
    features = [feat for feat in polygon_layer.getFeatures()]
    ids = []
    for buffers in buffer_layer.getFeatures():
        for polygons in polygon_layer.getFeatures():
            if buffers.geometry().intersects(polygons.geometry()):
                ids.append(polygons)
    intersect_layer = QgsVectorLayer("Polygon?crs=epsg:27700", 'Intersect', "memory")
    intersect_layer.dataProvider().addAttributes(polygon_layer.fields())
    with edit(intersect_layer):
        for feat in ids:
            feature = QgsFeature()
            feature.setGeometry(feat.geometry())
            feature.setAttributes(feat.attributes())
            intersect_layer.addFeature(feature)
    QgsMapLayerRegistry.instance().addMapLayer(intersect_layer)

func(point_layer, polygon_layer, 200)


texts: ArcMap Text Annotations Display in Proportion to Map Scale

texts:
My question is probably of an 'easy fix' however I've been struggling with it for a while! I'm trying to set up an annotation labels for several feature classess in a map and would like them to display proportionally - relatively to map scale. I use ArcGIS v 10.1 and Maplex label engine.

So - as an example, I am working on a regional map at a scale of 1:700K. My annotation on a FC (within GDB) is set to 12pt size. When I zoom in to a larger scale of 1:350K my annotations appear exactly the same size and are consequently tooo big!

I understand one can work with map document annotation or geodatabse annotation. I've seen some posts about setting a data frame reference scale to zero (eg http://forums.arcgis.com/threads/74004-Annotation-labels-at-different-scales), which should make mxd annotations display proportionally to map scale. However cannot seem to add custom reference scale of '0' to the customized list.

I'd greatly appreciate your comments on best practices or any other advices you may have!

Many thanks,

Magda


texts: Print Composer - raster layer name as text field

texts: I have three Landsat scenes and one polygon shapefile (column "region"). The atlas option is great to create prints of each "region". Is it possible to show the Landsat (raster) layer name automatically as a text field?


texts: Python problem during CADInput installation?

texts:
I was trying to install CADInput (version 0.60) on 2.4 Chugiak and 2.6 Brighton as well. In both cases I got an error message during installation. For me the following extract of that message might help someone to give me a hint to fix the problem:

UnicodeDecodeError: 'ascii' codec can't decode byte 0xdf in position 27: ordinal not in range(128)

Python-Version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

For some reason the installed version of Python does not work properly with the version of CADInput.

I need help to get it run. I'm not a programmer so maybe there's an easy way.


texts: How do you publish map caches from ArcGIS Data and Maps for Server?

texts:
I haven't published ESRI-supplied data before, but I followed the instructions that came with the Data and Maps for ArcGIS Server 10. When I view the supplied MXD file for a map cache, it's just a black background. The help docs say it's supposed to be different from what you expect, so I followed the rest of the steps:

Copy map cache folders to the arcgisserver\arcgiscache folder.
Copy the map documents to a shared location accessible to your GIS server. (Also copy the data folder from the documents folder).
Create a new map service using the MXD as the source document.

I followed the steps that come with the data but the service never displays any data. I suspect that the imprecise step of "Copy the map documents to a shared location accessible to your GIS server" requires something I haven't done. Has anyone published map services from the ArcGIS data discs before?


texts:
I am trying to create contours from a DEM in QGIS 3.4. I used the GDAL raster extraction contour tool. The contours are drawn correctly, but they are mirrored along the north line of the DEM (see image). Is there a setting that I missed when generating the contours? Both the CRS of the DEM and the CRS of the project are UTM NAD83 Zone 10 North.

Here is the output log:

Processing algorithm…
Algorithm 'Contour' starting…
Input parameters:
{ 'BAND' : 1, 'CREATE_3D' : False, 'FIELD_NAME' : 'ELEV', 'IGNORE_NODATA' : False, 'INPUT' : 'C:/Users/talldrin/Documents/Northstar/18-057 Cal Park Lakes/Support Docs/Lidar/Merged_DEM.ovr', 'INTERVAL' : 1, 'NODATA' : None, 'OFFSET' : 0, 'OPTIONS' : '', 'OUTPUT' : 'C:/Users/talldrin/AppData/Local/Temp/processing_a23c4b552e994e69b31a00e2657e2022/378e907d22c44a6f9c1389afc43c00ec/OUTPUT.shp' }

GDAL command:
gdal_contour -b 1 -a ELEV -i 1.0 -f "ESRI Shapefile" "C:/Users/talldrin/Documents/Northstar/18-057 Cal Park Lakes/Support Docs/Lidar/Merged_DEM.ovr" C:/Users/talldrin/AppData/Local/Temp/processing_a23c4b552e994e69b31a00e2657e2022/378e907d22c44a6f9c1389afc43c00ec/OUTPUT.shp
GDAL command output:
0...10...20...30...40...50...60...70...80...90...100 - done.

Execution completed in 1.35 seconds
Results:
{'OUTPUT': }

Loading resulting layers
Algorithm 'Contour' finished


texts: Add a legend in an ArcGIS Javascript API

texts:
I am quite new to Javascript. I have made a html script that includes services I published to ArcGIS for server. It all does exactly what I want and I made it by copying working scripts together and getting help from more experienced scripters. I want to finetune my service by adding a legend. Unfortunately I am not able to create a working script. I have added the parts of my script I assume that are necessary to create my legend. Does anybody see what I am doing wrong?

EDIT: I have edited the code. I am still not done. It does give me a legend on my browser but it says there is no content in there. On my arcgis server manager there is a legend.


texts: Google Earth Engine ROI area falling outside/partial coverage by Sentinel 1 Tile

texts:
Using this code I am trying to extract surface area of water bodies but some ROIs are not completely covered because of missing adjacent tile (see image). What can I do to add and mosaic the adjacent Sentinel 1 tile?
code link: https://code.earthengine.google.com/24e31dff3a145f38bbfd7deb59636da5
// Load region defined by polygon and add it to the map

Map.addLayer(roi, {}, 'ROI')
Map.centerObject(roi, 12)

//Load Sentinel-1 SAR collection and filter according to data collection type
var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(roi)
  .filterDate('2014-01-01','2020-07-29')
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))

//Add first image to map to get an idea of what a SAR image looks like
Map.addLayer(S1.first(),{bands: 'VV',min: -18, max: 0}, 'SAR image')

// Filter speckle noise
var filterSpeckles = function(img) {
  var vv = img.select('VV') //select the VV polarization band
  var vv_smoothed = vv.focal_median(100,'circle','meters').rename('VV_Filtered') //Apply a focal median filter
  return img.addBands(vv_smoothed) // Add filtered VV band to original image
}

// Map speckle noise filter across collection. Result is same collection, with smoothed VV band added to each image
S1 = S1.map(filterSpeckles)

//Add speckle filtered image to map to sompare with raw SAR image
Map.addLayer(S1.first(),{bands: 'VV_Filtered',min: -18, max: 0}, 'Filtered SAR image')

//Classify water pixels using a set threshhold
//Here we are using -16. This is only an approximation and will result in some errors. Try adjusting the
var classifyWater = function(img) {
  var vv = img.select('VV_Filtered')
  var water = vv.lt(-16).rename('Water')  //Identify all pixels below threshold and set them equal to 1. All other pixels set to 0
  water = water.updateMask(water) //Remove all pixels equal to 0
  return img.addBands(water)  //Return image with added classified water band
}

//Map classification across sentinel-1 collection and print to console to inspect
S1 = S1.map(classifyWater)
print(S1)

//Make time series of water pixels within region
var ClassChart = ui.Chart.image.series({
  imageCollection: S1.select('Water'),
  region: roi,
  reducer: ee.Reducer.sum(),
  scale: 100,
})
  .setOptions({
      title: 'Inundated Pixels',
      hAxis: {'title': 'Date'},
      vAxis: {'title': 'Number of Inundated Pixels'},
      lineWidth: 2
    })

//Set the postion of the chart and add it to the map
ClassChart.style().set({
    position: 'bottom-right',
    width: '500px',
    height: '300px'
  });

Map.add(ClassChart)

// Create a label on the map.
var label = ui.Label('Click a point on the chart to show the image for that date.');
Map.add(label);

//Create callbakc function that adds image to the map coresponding with clicked data point on chart
ClassChart.onClick(function(xValue, yValue, seriesName) {
    if (!xValue) return;  // Selection was cleared.

    // Show the image for the clicked date.
    var equalDate = ee.Filter.equals('system:time_start', xValue);
    //Find image coresponding with clicked data and clip water classification to roi
    var classification = ee.Image(S1.filter(equalDate).first()).clip(roi).select('Water');
    var SARimage = ee.Image(S1.filter(equalDate).first());
    //Make map layer based on SAR image, reset the map layers, and add this new layer
    var S1Layer = ui.Map.Layer(SARimage, {
      bands: ['VV'],
      max: 0,
      min: -20
    });
    Map.layers().reset([S1Layer]);
    var visParams = {
      min: 0,
      max: 1,
      palette: ['#FFFFFF','#0000FF']
    }
    //Add water classification on top of SAR image
    Map.addLayer(classification,visParams,'Water')

    // Show a label with the date on the map.
    label.setValue((new Date(xValue)).toUTCString());
  });


texts: Identifying "chains" of adjacent polygons with PostGIS

texts:
I have a table of polygons (land parcels in a city) and I am trying to identify "chains" of 5 or more adjacent polygons. For example, if each letter is a polygon and each . is empty space, AB.CD.EFG.HIJKL.MNOP.QRSTU.VWXYZ then I want HIJKL, QRSTU and VWXYZ. In real life they aren't in a line like that of course.

The messy solution I used for a first attempt is this query:

SELECT p1.parcel_c as "Parcel 1", p2.parcel_c as "Parcel 2", p3.parcel_c as "Parcel 3", p4.parcel_c as "Parcel 4", p5.parcel_c as "Parcel 5", ST_Union(array[p1.geom, p2.geom, p3.geom, p4.geom, p5.geom]) as merged_geom FROM abandonded_vacant_lots_table p1 INNER JOIN abandonded_vacant_lots_table as p2 ON ST_TOUCHES(p1.geom, p2.geom) INNER JOIN abandonded_vacant_lots_table as p3 ON ST_TOUCHES(p2.geom, p3.geom) INNER JOIN abandonded_vacant_lots_table as p4 ON ST_TOUCHES(p3.geom, p4.geom) INNER JOIN abandonded_vacant_lots_table as p5 ON ST_TOUCHES(p4.geom, p5.geom) WHERE p1.parcel_c not in(p2.parcel_c, p3.parcel_c, p4.parcel_c, p5.parcel_c) and p2.parcel_c not in(p1.parcel_c, p3.parcel_c, p4.parcel_c, p5.parcel_c) and p3.parcel_c not in(p1.parcel_c, p2.parcel_c, p4.parcel_c, p5.parcel_c) and p4.parcel_c not in(p1.parcel_c, p2.parcel_c, p3.parcel_c, p5.parcel_c) and p5.parcel_c not in(p1.parcel_c, p2.parcel_c, p3.parcel_c, p4.parcel_c)

As you can see this is clumsy and inelegant and I suspect inefficient. It also returns a lot of almost duplicates when there is a chain of more than 5 parcels, eg if the parcels are ABCDEF it returns ABCDE and BCDEF when a better solution would return the longest chain, ABCDEF. How can I make this query more efficient, more elegant and in general more versatile?


texts: Mistake in Calculating Area or a Bug in QGIS?

texts:
I am trying to calculate the area of a polygon following the step by step mentioned in all the questions related to this issues in this website. Nonetheless, the result of the area calculation is about 6.1 and of course it is wrong cause I was expecting something more than 50000 km². I don't know if I am mistaking in some point or there is a bug in the qgis version I have.
I have:

changed the layer projection to UTM (In this case is UTM 23);
changed the project layer projection to UTM;
calculated in attributes/field calculator/geometry/$area.

Still the result is something weird. This is the polygon I am working Três Marias Polygon


texts: Creating polygon of difference using ArcMap

texts:
I have two polygon layers, each describing the Geomorphic units (GU) within a river channel across two surveys.
Using ArcMap (10.7.1), I would like to create a new layer, illustrating the change in the number and area of the different geomorphic units between the two surveys. I have attributes describing area, GU's etc. The number, classification and shape of the polygons differ between surveys, they do not share FID's.

I would like to create a .csv or (ideally) a shapefile showing

which GU's changed and what they changed from and into (i.e. pool to riffle).
Union tool does not produce the outcome I want; it creates a new shapefile of both surveys merged together.


texts:
I trying to develop a PHP script that will simply applyedits to a feature services.

Right now I'm trying to use the In the ESRI restend point example. I'm thinking if someone can help me, point me to some doc, show me how to use the esri endpoint.  I will be able to take that script learn from it and move from there.  There are several moving parts I'm hoping to accomplish but right now I'm stuck on this piece.

If I hit this esri rest endpoint:

  sampleserver3.arcgisonline.co...r/0/applyEdits

and add these parameters to the add input box

  [{"geometry":{"x":-13085962.2699721,"y":4061306.82523268,"spatialReference":{"wkid":102100}},"attributes":{"description":"Test Nat","type":5}},{"geometry":{"x":-13072136.6421672,"y":4073467.58686212,"spatialReference":{"wkid":102100}},"attributes":{"description":"klsf","type":0}}]

I can add a feature.

I want to create a simple PHP script that will do this for me.  CAN someone show how to do this?

If I can get this working I'm hoping I can feed the parameters from an xml feed and consistently push data to my rest points. Thanks.


texts: Calculating average pixel value for each band across multiple raster layers?

texts:
I have a set of raster files that are essentially data sets of precipitation across a region. I have multiple raster files that are the outputs of different climate models. Each pixel in the raster file has 365 bands that represent the precipitation value for each day. I want to calculate the mean for each pixel for each band across the different raster files.

I tried converting the raster to point files and extracting multiple values to the point, extracting the data to excel and then averaging them for each raster. I was wondering if there is a more efficient way to do this.


texts: Make IDataElement from in memory featureclass

texts:
I'm trying to do some fieldmapping, but need to make a dataelement from an in-memory featureclass in ArcObjects to pass to the geoprocessor.

Basically, this line from the sample code:

IDETable inputTableA = (IDETable)gputilities.MakeDataElement(@"C:\Temp\k\tmpExplode.gdb\ExplodeOUTPUT", null, null);

But I there is no "real" path to the featureclass.  Also tried:

IDETable inputTableA = (IDETable)gputilities.MakeDataElementFromNameObject(null);

but I can not get to an IName for it either.  Has anyone made an IGPFieldMapping object from an in-memory featureclass before.  My code works fine with a file geodatabase layer.


texts: How to stylize/design maps without tools like m*p b*x (paid tools) using OSM data

texts: I would like to know, how I'm able to stylize/design maps without using paid tools. I would love to work with Openstreetmap data. I know I could stylize these data in free GIS software like QGIS but after uploading to server I would like to add some additional functionalities using JavaScript. Any suggestions?


texts: Batch Mosaic Rasters: arcpy.MosaicToNewRaster_management (ERROR 000628: Cannot set input into parameter coordinate_system_for_the_raster)

texts:
I am trying to batch mosaic rasters using arcpy.MosaicToNewRaster_management in arcpy:

    #This is for mosaicking a lot of rasters
    import arcpy
    from arcpy import env
    arcpy.CheckOutExtension("Spatial")
    env.overwriteOutput = True

    #Set the current workspace
    env.workspace = r"C:\thesis\for_sampling\sampling_outputs\new_loss"
    env.nodata = "MINIMUM"
    env.compression = "LZ77"

    coordinate_system =  'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],\
                       #UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG",4326]]'

    list_all_rasters = ["newlossGreater_Luzonforests_onlyloss_20N_120E.tif;newlossGreater_Luzonforests_onlyloss_30N_120E.tif;\
newlossGreater_Luzonforests_onlyloss_20N_110E.tif","newlossGreater_Palawanforests_onlyloss_20N_120E.tif;\
newlossGreater_Palawanforests_onlyloss_10N_110E.tif;\newlossGreater_Palawanforests_onlyloss_20N_110E.tif",\
"newlossGreater_Negros_Panayforests_onlyloss_10N_120E.tif;newlossGreater_Negros_Panayforests_onlyloss_10N_110E.tif",\
"newlossGreater_Mindanaoforests_onlyloss_10N_120E.tif;newlossGreater_Mindanaoforests_onlyloss_10N_110E.tif;\
newlossGreater_Mindanaoforests_onlyloss_20N_110E.tif"]

    output_list = ["Luzon_loss.tif","Palawan_losss.tif","Negros_Panay_loss.tif","Mindanao_loss.tif"]

    for raster in  range(0,4):
       output_list = ["Luzon_loss.tif","Palawan_losss.tif","Negros_Panay_loss.tif","Mindanao_loss.tif"]
       env.workspace = r"C:\thesis\for_sampling\sampling_outputs\new_loss"
       print raster # checking the list
       arcpy.MosaicToNewRaster_management(list_all_rasters[raster], r"C:\thesis\for_sampling\sampling_outputs" , output_list[raster], 1, coordinate_system)
       print raster+ " is ok!"
    print "Finish all of them!"

The error talks about coordinate systems. I've tried to use "GCS_WGS_1984.prj" but it didn't work. I also put that .prj file into the current workspace, still the same error:

    arcpy.MosaicToNewRaster_management(list_all_rasters[raster], r"C:\thesis\for_sampling\sampling_outputs" , output_list[raster], 1, coordinate_system)
    File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 11613, in MosaicToNewRaster
raise e
    ExecuteError: ERROR 000622: Failed to execute (Mosaic To New Raster). Parameters are not valid.
    ERROR 000628: Cannot set input into parameter coordinate_system_for_the_raster.

I also tried to change the variable: coordinate_system = arcpy.SpatialReference(4326), but the error persists.


texts: Trying to import Tiff files created in arc9.3 into 10.2, not working

texts:
I'm trying to use the scanned maps available here: http://www.env.gov.bc.ca/esd/distdata/ecosystems/TEI_Scanned_Maps/B01/B01-4700/

I save all 6 files in a folder, open ArcGIS 10.2 and try to add to my data frame from the catalogue. I can view the tiff but it will not georeference. Warning shows "Unknown Spatial Reference". This should be embedded in the provided files but are not being accessed properly.

Can someone try to import the files from the link above and let me know if they can get them to properly georeference?


texts: Why UPDATE from Analysis tools in ArcGIS 10.3.1 is not available to me, does it need some specific license or can it be download separately? How is it possible to achieve exactly the same result as with UPDATE without it?


texts:
I am looking for a tool or query to split the OSM highways into multiple segments. Each segment should only consist of two nodes (start and end).

I should do the same as described in this GIS Stackexchange answer but with PostGIS.

Is there an importing tool like osm2po that I can do this with?
Note I don't only want to split it at intersections.

Or alternatively, how could I use ST_Split() or ST_Dump() to do this?


texts: Why does projection WKT differ between osr.SpatialReference and gdalinfo

texts:
I am writing a simple geotiff using the GDAL python bindings, but the projection WKT reported by gdalinfo differs from the source. The script

from osgeo import gdal, osr
import numpy

# Create synthetic data
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
raster[256, 256] = 1
top_left = (100000, 200000)
pixel_size = 30

output_fn = 'test.tif'
driver = gdal.GetDriverByName('GTiff')
geotiff = driver.Create(output_fn, 512, 512, 1, gdal.GDT_Byte)
src_geotran = [top_left[0], pixel_size, 0,
               top_left[1], 0, -pixel_size]
geotiff.SetGeoTransform(src_geotran)

srs = osr.SpatialReference()
srs.ImportFromEPSG(3413)
srs_wkt = srs.ExportToPrettyWkt()
print(srs_wkt)

geotiff.SetProjection(srs_wkt)
geotiff.GetRasterBand(1).WriteArray(raster)
geotiff = None

produces the WKT

PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",70],
    PARAMETER["central_meridian",-45],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","3413"]]

Whereas the WKT from gdalinfo is

PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
    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["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3413"]]

Notably, the latitude_of_origin and central_meridian show as 0 in gdalinfo, when they should be 70 and -45 respectively.

Tested with Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1


texts: Ubuntu QGIS 2.4.0 SSL connection to remote PostGIS fails b/c can't specify certificates

texts:
I can connect with SSL in PGAdmin to my remote server, but it requires specifying the filepaths for certificate files due to my having different certs for different servers.

However, when connecting from QGIS, I cannot find the option to specify the certificates to use to connect to my remote host. So naturally, when I try connecting and require SSL I get

  SSL error: certificate verify failed


texts: Pixelated result using SVG as label background?

texts:
I have points that are labelled with a 2-digit numeric attribute, shown by default with a circle background (white fill, coloured outline of 1 map unit, buffer of 1 map unit beyond the label).

I have been asked to show deleted points (deleted since the last iteration of the map) with a diagonal slash through the number. I am attempting to do this using an SVG as the label background, but the results are pixelated. The width of the stroke element of the SVG also seems to change with zoom level. I am not offered an option for the units of stroke width on the SVG background, just "symbol units".

Are there settings I should check either in the original SVG file or in QGIS that might help with either or both of these?

I am using this SVG: https://commons.wikimedia.org/wiki/File:Cross_within_circle_2.svg
- I do not know why only one of the diagonal lines shows up (one or both would be fine). I don't know SVGs well but I am familiar with markup languages in principle and happy to attempt edits if provided with a little guidance.

Alternatively, is there another way that I can have the labels mimic the non-deleted points, just with the addition of a diagonal line?

I have suggested just using strikethrough text but this is deemed not clear enough by the recipient. However, at present it looks like a better alternative if I cannot fix the SVG issue.

(QGIS 3.4.7, Win10)


texts: Determining lake catchment in Arc Hydro tools?

texts:
I am trying to determine a lake's catchment with Arc Hydro tools or other tools in ArcMap.

I tried basin tools in sptial analyst and catchment grid delineation tools in ArcHydro Tools but all the time I got too large areas along the watercourses where there are many lakes. I need the catchment of individual lakes. How to do it?


texts: How to change the styles (fonts, in particular) in a ArcToolbox tool UI?

texts:
I am attempting to change the font size in my ArcToolbox tool UI. I have scoured the ArcGIS 10.0 online help and found nothing on this topic. Googling it turned up the 9.1 version, which has been helpful. I followed those instructions, and here is what I have done:

I copied MdDlgContent.xsl that styles Tools to a local directory, the same as my Tool actually
I edited the CSS and the HTML in MdDlgContent.xsl and have been successful in making a quite ugly Tool (just for the purpose of seeing whether I can make changes occur or not) where the background color, table border, table cell padding, and background image are changing:

Now, that said, I cannot, for the life of me, get the caption fonts (surrounded in red boxes above) to change. I have added/edited CSS and HTML all over the code, trying every dirty non-validating hack I know of. I would post the code, but it's just too much. The XSL, on my system, is located at:

C:\Program Files (x86)\ArcGIS\Desktop10.0\ArcToolbox\Stylesheets

If you scroll down in the 9.1 help on this topic, they show an example where the font size and color have both been changed.

I'm thinking this is something that is right in front of me that I'm just not seeing and am wondering if anyone has figured this out before.


texts:
I need a detailed description of a GPS trajectories. I don't need to know what a GPS trajectory is, I already know it. Exactly I need a paper (better if scientific cut) that talk about different type of trajectories. e.g.: Difference in pedestrian, marine and car  and highlight what difference presents from each other, can anyone help?

Best regards


texts: Defining polygon's area from OpenStreetMap

texts:
In OpenStreetMap, I have a multipolygon OSM feature with its relation ID 6040040 (example).

How do I calculate the area of a such features using QGIS?


texts: Why can't I edit layers created with DBManager?

texts:
I noticed a strange behavior in layers created with DB manager.
Tested it in QGIS 2.6.1 and QGIS 2.12.0

In QGIS 2.6.1.
When i create a layer based on a sql query with the setting "Avoid selecting by feature id" turned off, i cannot start editing session on this new layer. The "start editing" button is inactive.
When i create a layer based on the same sql query with "Avoid selecting by feature id" turned on, the editing button in QGIS gets activated, i can start editing session but "Add feature" and "Save layer edits" buttons stay inactive. I can move any features or nodes on a feature, but still i cannot save the edits.

In QGIS 2.12.0
There is no difference, if i set "Avoid selecting by feature id" on or off, editing button is inactive and i cannot start editing session.

Is it a bug, am i doing something wrong or is QGIS ment to behave like this?
And what is this "Avoid selecting by feature id" option for?

O, and the data i want to edit comes from a PostGis database.


texts:
I have several installs of Postgres / PostGIS, version 9.4 up to and including 9.6. I would like to upgrade all the databases (maybe 20 total) in the various installs to a single install of the current version.

How do I migrate the databases to the current version? Do I just backup each database in situ, restore in the new PG instance and then upgrade the PostGIS extension (like ALTER EXTENSION postgis UPDATE TO "2.3.1";)? Any other steps?


texts: ArcGIS Server 10.0 too slow to generate cache tiles

texts:
I want to define the zoom level scale for my arcgis web map, but it seems that I should build the cache tiles first (in ArcMap). I have waited for 2 hours but the cache status is still maintain in 0%. I was wondering if this is normal? Thanks.

Details:

All files in TIFF (elevation map): (28,4 GB)
ArcGIS Server 10.0 Enterprise Windows 7 32 bit.
In Analyze Map 10 Warnings and 7 Messages
7 scales

How to speed it up?


texts:
I'm using openlayers with OSM as a map source. I would like to change map settings co I cannot zoom more than map size.

This is a situation now. I do not want to see the white space. I do not want map to wrap X (I already did that with wrapX : false in OSM constructor). But how can I set up the map to see it just once? I tried to set min zoom to some value But I could still see the white space.

This is my init function

    this._source = new OSM({
        cacheSize: 10000,
        wrapX : false
    });

    this._layer = new OlTileLayer({
        source: this._source
    });

    this._view = new OlView({
        center: [-7916041.528716288, 5228379.045749711],
        zoom: 4,
        minZoom: 2
    });

    this.map = new OlMap({
        target: 'map',
        layers: [this._layer],
        view: this._view,
        controls: [
            new OlControlZoom(),
            new OlControlZoomSlider()
        ]
    });

Can anyone help me?


texts: QGIS labeling: expression that shows values from one attribute except when the objects has a value x at another attribute

texts:
I want to label each object of my layer with the values of a certain attriubte. However, if an object has the value "1" at another attribute, this label should be invisible or not rendered.

Any ideas how to do that?


texts: Polyline buffer based on attribute value in ArcMap, choice in join/cap style?

texts:
In QGIS I had no difficulty to create a buffer (in fact a single symbol) whose width is based on an attribute value of the polyline vector:

In ArcMap, I don't find any good solution. I tried:

the buffer tool (in Analysis Tools/Proximity/Buffer), which doesn't
let me choose the symbol and the side/end type,
and the symbology > Multiple attributes > Symbol Size with the relevant attribute and "template" symbol, which doesn't apply the style chosen. The result is instead a green polyline whose width varies according to the attribute, but the style is not applied, as shown below (see the green line in the background, even if I had chosen a black/purple template):

Am I wrong somewhere?


texts: Converting between Esri Geometry and WKT using ArcObjects?

texts:
It seems like there is no means of converting an ArcObjects geometry to the Well-Known Text representation (and vice versa) in ArcGIS ArcObjects API. The only thing I was able to find is conversion to WKB (the IWkb interface).

Is there a way to perform conversion between geometry objects and WKT or do I have to implement it myself? I primarily aim for .NET implementation without much external dependencies.


texts: How to calculate the azimuth of a line of a multilinestring in PostGIS?

texts:
I have two spatially enabled tables in PostgreSQL with a geometry in each:
First table is with POINT type, second is with MULTILINESTRING.

I need to calculate the azimuth of the LINESTRING part of a MULTILINESTRING, which is closest to the POINT.

So the question is: how do I extract the LINESTRING in question, and compute the azimuth?


texts: midpoint of segments in QGIS determination

texts:
This is a follow up to How to find the middle point of a line in QGIS - Where to place input file name

Errors coming up when I executed the code below (taken from Finding middle point (midpoint) of line in QGIS?) included an "IndentationError: unexpected indent" which I resolved by deleting a leading space on the 8th line of the code right before for i in range(1, len(list)):.

A second set of errors came up "Traceback (most recent call last): File "", line 1, in  NameError: name 'mylayer' is not defined". I think I resolved the 'mylayer' error part by providing the name of the active layer ('col19_lines') on line 21 of the code for elem in 'mylayer'.selectedFeatures(): I want midpoints computed from.

The correction I made still did not give me the midpoints as another error came up "Traceback (most recent call last): File "", line 1, in  AttributeError: 'str' object has no attribute 'selectedFeatures'".

How can I make Gene's code run without error?

QGIS 1.8 is being used.

code:

def mid(pt1, pt2):
   x = (pt1.x() + pt2.x())/2
   y = (pt1.y() + pt2.y())/2
   return QgsPoint(x,y)

def pair(list):
   '''Iterate over pairs in a list '''
   for i in range(1, len(list)):
        yield list[i-1], list[i]

def create_geometry(point,pr):
    # create geometry record
    seg = QgsFeature()
    seg.setGeometry(QgsGeometry.fromPoint(point))
    pr.addFeatures( [seg] )

# memory layer
pt_layer = QgsVectorLayer("Point", "midpoint", "memory")
pr = pt_layer.dataProvider()

for elem in 'col19_lines'.selectedFeatures():
    line = elem.geometry()
    for seg_start, seg_end in pair(line.asPolyline()):
       line_start = QgsPoint(seg_start)
       line_end = QgsPoint(seg_end)
       # midpoint
       midpt= mid(line_start, line_end)
       # add midpoint point to layer
       create_geometry(midpt,pr)
       pt_layer.updateExtents()

QgsMapLayerRegistry.instance().addMapLayers([pt_layer])

I need to properly process a layer named 'col19_lines' I do not have a background in python scripting. What do I need to do for the routine to process specifically a layer named 'col19_lines' so the midpoint of segments will be computed and no "Traceback (most recent call last): File ", line 1, in AttributeError: 'str' object has no attribute 'selectedFeatures'" error shows up.


texts:
I have data of monthly temperature for 60 years.

I want to analyze the trend for the monthly temperature in R.

How do I construct R Code which would analyze monthly trend over the time series to determine spatial distribution of significant trends in monthly average and min/max data? I have this code but do not know how to adjust so

//get the direction and magnitude of trends,
time <- 1:nlayers(gimms.sum)
fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[2] }}
gimms.slope=calc(gimms.sum, fun)
gimms.slope=gimms.slope*31 //number of years
plot(gimms.slope)
fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[8] }}
p <- calc(gimms.sum, fun=fun)
plot(p, main="p-Value")
//mask all values >0.05 to get a confidence level of 95%:
m = c(0, 0.05, 1, 0.05, 1, 0)
rclmat = matrix(m, ncol=3, byrow=TRUE)
p.mask = reclassify(p, rclmat)
fun=function(x) { x[x<1] <- NA; return(x)}
p.mask.NA = calc(p.mask, fun)
//Mask insignificant trends
trend.sig = mask(gimms.slope, p.mask.NA)
plot(trend.sig, main="significant NDVI change")


texts:
I am trying to make a map depicting forests in Japan using TIF files downloaded here. I want to reproject the data, clip it, and use it to style a hillshade map. When I attempt to reproject it however, it takes a very long time, and the output file is huge! The original TIF is quite big at 411.6MB, however the reprojected file clocks in at 1.75GB, and this is only one of roughly 8 TIF files I eventually want to use.

Why is this happening? Is there some fundamental step that I'm missing? If there's any other information I can provide, please let me know.

Data for the original file:

Name    1_Hansen_GFC-2018-v1.6_treecover2000_50N_130E
CRS EPSG:4326 - WGS 84 - Geographic
Extent  130.0000000000000000,40.0000000000000000 : 140.0000000000000000,50.0000000000000000
Unit    degrees
Width   40000
Height  40000
Data type   Byte - Eight bit unsigned integer
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description /Users/benleamon/Documents/Infographics/projects/japan-forests/data/forests/1_Hansen_GFC-2018-v1.6_treecover2000_50N_130E.tif
Compression LZW
Band 1
More information
AREA_OR_POINT=Area
Dimensions  X: 40000 Y: 40000 Bands: 1
Origin  130,50
Pixel Size  0.0002500000000000000052,-0.0002500000000000000052

And then here is the data for the output file:

Name    reproj-full-forests
CRS EPSG:6691 - JGD2011 / UTM zone 54N - Projected
Extent  -439944.6735999999800697,4428235.8196999998763204 : 428345.5341999999945983,5596865.5570999998599291
Unit    meters
Width   36028
Height  48490
Data type   Byte - Eight bit unsigned integer
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description /Users/benleamon/Documents/Infographics/projects/japan-forests/output/reproj-full-forests.tif
Compression
Band 1
More information
AREA_OR_POINT=Area
Dimensions  X: 36028 Y: 48490 Bands: 1
Origin  -439945,5.59687e+06
Pixel Size  24.10042766181858553,-24.10042766343575948

My system is a macbook air with 4GB ram, which I know is far less than ideal, hence my plan to process and clip each tile independently before combining on the final map.


texts: Access to ArcGIS Online data without REST API

texts: is there any possibility to get information/values from hosted ArcGIS Online tables for a locally running python script without using the ArcGIS REST API?


texts: Out of stack - Internet Explorer - WebAppBuilder

texts:
I'm developing a viewer with ESRI WebAppBuilder. I have a widget that use a table programatically created with component dojox.grid.DataGrid.
When Internet Explorer (version 11) renders datagrid with this code grid.startup(), I am getting "Stack overflow at line 0".
Do you know how to debug/fix this?


texts: Sentinel-2 imagery granules downloaded from AWS archives (http://sentinel-pds.s3-website.eu-central-1.amazonaws.com) have a little different data structure than the original granule sub-folders of ESA SAFE product. Consequently, it seems the SNAP tool seems not able to import the data as one product with metadata - and I need to process the data to level 2 (atmospherically corrected ground reflectance). Is there a way? Or is there a site which allows downloading individual granules in original SAFE data structure?


texts: How does TomTom store trackpoint information on its internal drive?

texts: How does TomTom stores tracks and trackpoints on its internal drive? I maintain an open source computer forensics tool and would like to extract the data from both files and unallocated space. My experience with Garmin GPS devices is that they all use the same internal data format, and I suspect that TomTom devices do the same.


texts: PostgreSQL - attach to point table the maxium distance value to edge of containing polygon

texts:
I'm trying to find the furthest distance from a point to the edge of the polygon it sits in. The polygons in this case have been converted to vertex points.

I have two POINT tables in Postgres (and a POLYGON table with 3 shapes).

One POINT table (vertextable) contains vertex points of the 3 polygons (about 18 points in total). I created this vertex point table using ST_DumpPoints.

The other POINT table (insidepointtable) contains 3 points, each one lies somewhere within each of the original polygons. Each point shares a unique ID number with the polygon vertex points they fall inside (see pic).

I want to create a new table containing the 3 points from insidepointtable and the maximum distance value of the point from vertextable.

I'm having trouble returning just the 3 distinct rows. When I use:

CREATE TABLE public.vertexmaxdistances AS
SELECT a.id, ST_MaxDistance(a.geom, b.geom)
FROM public.insidepointtable a, public.vertextable b
WHERE a.id = b.id
ORDER BY id DESC, st_maxdistance DESC;

I get every distance to every point returned. How can I get:

id   |  st_maxdistance
______________________
1    |             5.5
2    |             4.2
3    |             4.9

(for example)

Furthermore, I'd actually like to do this distance attachment for millions of sites. So ideally the solution needs to be scaleable with decent performance.


texts:
I'm just a rookie when it comes to creating models in ArcGIS. I have two shapefiles from which I need to calculate an AVG and a MAX erosion rate for individual properties.

The two shapefiles are:

Coastal_Properties & Coastal_Change

Coastal_Properties has a Property ID (PID), and Coastal_Change has a Rate of erosion in metres per year (RChngMYr).

I'd like to retrieve an Average and a Max RChngMYr for each individual PID in the data-set.

I have no problem running this query using SQL select in MapInfo, but am unsure how to go about it in ArcGIS. It will be better in the long run to create a model in Arc as Property boundaries are dynamic, and it would be terrible to have to run the query in MapInfo everyday. I'd like to create a Geoprocess service from this model and use it in a web application where it would return a phrase such as "The Max erosion rate for PID 12345 is 3.1 m/yr".

Any help is appreciated! Thanks in advance


texts: Saving project creates three layers per one original layer and some files are "no valid data sources"

texts:
It doesn't matter how I save my project (Qconsolidate-Plugin, Relocator-Plugin, QPackage-Plugin or manual), the output is always defective: The export always consists (as shown in the picture) of three layers per one original layer in my original project: One .shp-layer, one .shx-layer and one .dbf-layer.
Additionally the .qgj-file, the .prj-file and the .cpg-file of every layer are "no valid data sources" referring to the defect report of QGIS.
The result is a wrong projection and a wrong symbology.

How ist that possible?


texts: Extracting polygons that overlay other polygons in different layer

texts: I have two polygon layers that have multiple polygons in them. One layer is an administrative boundary layer and the other is a layer showing crop cover types. I would like to extract only administrative boundaries that have the certain crop type within them. Is there a way to do this in QGIS?


texts: Authentication with ArcGIS Server Identity

texts:
Just getting back into ArcGIS for .NET again after about 2 years of doing the training course.

I have a question about the "Add ArcGIS Identity" that you choose in Visual Studio - it doesn't seem to matter what combination of username/password I think of, it always says I have an invalid username/password (by the way I am using accounts that are in the agsadmin group) on the AGS server. Yet, if I copy a ciphered identity tag from the web.config of some random AGS application I made via ArcGIS Manager - it works, i.e. the web app loads up the map services and I can see a map in my application.

Now, I can't decipher this identity tag as you can probably guess - but does anyone know what I might be doing wrong or what I need to look, or check for in order to get the correct credentials right. I just can't rest until I have this figured out - instead of doing this hack.

UPDATE - Answer

I made sure my dev PC had the ArcGISWebServices user set with the same password as the AGS physical server - so the "impersonation" was consistent between the two machines and accounts. I replaced out the ciphered details in the web.config and replaced with the ArcGISWebServices one and it works fine. I also tried doing it via the "Add ArcGIS Identity" and now it works but make sure you put the machine name (i.e. the name of your dev pc) in the "Domain or Machine Name" textbox.


texts: Convert OMI/AURA No2 data from molecules/cm2 to micrograms/m3 (µg/m3)

texts: Currently, I am working on OMI/AURA NO2 data which is in molecules/cm2. I need to convert it to micrograms/m3 (µg/m3) to compare with the O3 concentrations. I have gone through the user guide for OMI data, but couldn't find anything related to the units & conversion.


texts: Test if a feature is visible in QGIS expression?

texts:
I'd like to test whether a given feature is displayed, i.e. visible on the map canvas, inside a QGIS expression. Any way to do this?
I'm essentially looking for an analog to the is_selected() function, which tests whether a feature has been selected, but querying whether it's displayed or not. For my use case, I'm more concerned about whether it meets the user's current filtering criteria in the layer tree (e.g. which categories are set to display in the categorized renderer) than whether it is within the displayed canvas extent, but I could live with both. I don't see such an analog in the list of functions in the QGIS documentation, however.
As further context, I'm building an expression that aggregates across points inside polygons. Simplified version is
aggregate(layer:='YetiSites',
aggregate:='sum',
expression:=NumberOfYetiObserved,
filter:=intersects($geometry,geometry(@parent)))

which adds up the number of yeti observed at all sites that within the zone defined by the current (polygon) feature. See below, where it's assumed this expression is used in the labels for the zone polygons.

However, the YetiSites layer also has an attribute Observer, and in the layer tree users can choose which observers' sightings they want to see on the map. If the display settings exclude, say, Bill's observations, I want my sum to exclude them too. So I'd like to be able to write something like
filter:=(intersects($geometry,geometry(@parent)) and is_displayed($currentfeature))

but I'm trying to find how to code the (ficticious) is_displayed(...) part. Without this clause, I get the following, with 6 rather than the desired 5 shown for Zone A since the aggregation picks up Bill's point even though it is not being displayed.

I can of course emulate this by running a processing algorithm, or defining the display filter upstream and then using it both when displaying YetiSites and in my aggregation expression. But I'd love to just be responsive to what the user is selecting in the layer tree to display.


texts: Can ArcGis use imagery from Google or Bing in order to compute Service Areas?

texts:
I am trying to make service areas and I am wondering if I can use one of the open source data for it, not as a basemap, but if the software can take the real-time information from these maps, such as oneway restrictions, speed limits and so on.

When working with normal layer for roads I don't have all this data and also it is not provided the history data so that I can use different time when computing the service areas.


texts: Creating records in output shapefile using Python

texts:
I have a list (output from a previous loop) containing attributes from an input shapefile. What I am stuck with is inputting those attributes to a new shapefile.

I have plotted the point data using w.point(x, y) and also created the relevant fields w.field() from the input shapefile fields.

Where I am stuck is with the copy of the records, which is giving me an Error related to Deletion Flag.

  AttributeError: 'tuple' object has no attribute 'startswith'

I am pasting a section of the code with my queries.

from datetime import datetime
import osgeo.ogr, osgeo.ogr
from osgeo import ogr
from osgeo import gdal
import shapefile
import os

sf = shapefile.Reader("-- Input File Location --") # Reading the Shapefile
fields = sf.fields # Reading the attribute fields
records = sf.records() # Reading the records of the features
shapRecs = sf.shapeRecords() # Read Geometry and Records simultaneously

w = shapefile.Writer(shapefile.POINT)
w.autoBalance = 1
result1 = []
result2 = [] # ??? Is the error because we cannot query for list in a list ???

print '** Call function pt(p) where p is node index... **' # Prompting the user to
                                                           # call the function with
                                                           # the required number of nodes
def pt(p):
    for i in range(len(shapRecs)):
        u = shapRecs[i].shape.points[p-1]   # Gives the XY Coordinates of
                                            # the queried points

        v = shapRecs[i].record[0:]          # Generating the Attributes
        result1.append(u)
        result2.append(v)

        w.point(u[0], u[1]) # Create points from coordinates each time the loop runs
        w.autoBalance = 1

    # -- Adding the fields
    w.field(fields[0]) # Adding the first field, WITH Deletion Flag (???)
    for i in range(1, len(fields)): # Loop begins from the second field (???)
        w.field(fields[i])
        w.autoBalance = 1

    # CODE SHOWS THE ERROR STATED ABOVE FROM THE w.record() statement
    # -- Adding the records
    for i in range(len(result2)):
        w.record(*result2[i])
        w.autoBalance = 1

    w.save('-- Output File Location --')


texts: QGIS standalone plugin: Display a progress bar with a busy state

texts:
I have a standalone QGIS plugin and to inidicate that the plugin is processing, I want to add a progress bar that indicates a busy state, but not showing the exact progress like in the picture below:

I use the following code:

def main(self):

    self.progress()

    ...some more code here...

    self.dlg.close()

def progress(self):

    self.dlg.progressBar.reset()
    self.dlg.progressBar.setRange(0,0)

So, the idea is the call the def progress at the beginning of my code and execute everything else afterwards and in the end the window closes. For the progress bar I used a minimum of 0, a maximum of 100 and a value of 0 in QT Designer, just to display an empty progress bar.

However, if I set the range to (0,0) (which should add a busy state to the progress bar), nothing happens. Any idea what I do wrong?


texts: Entering date with arcpy.GetParameterAsText()?

texts:
I want to allow users to enter a date as a parameter through a dialog window in a toolbox script.

I have obviously defined the field with a "Date" data type but something is wrong with the formatting.

arcpy.CalculateField_management(fc, "Date", "01.07.2015", "PYTHON", "")

If I enter "01.07.2015" directly in the FieldCalculator in ArcMap it works, but I can´t define it as a parameter in arcpy? What am I doing wrong?

import sys
import arcpy
from arcpy import env

# allow overwriteOutput
env.overwriteOutput = True

# for Script
FC = arcpy.GetParameterAsText(0)
Date = arcpy.GetParameterAsText(1)

try:

    arcpy.CalculateField_management(fc, "Stand", Date, "PYTHON", "")

except Exception, e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message

This is in Modelbuilder. The code runs, but the field is not updated. If I take remove the modelbuilder parameters and enter them in the python script the error is:

import sys
import arcpy
from arcpy import env

# allow overwriteOutput
env.overwriteOutput = True

Date = "01.01.2001"

try:

    arcpy.CalculateField_management(fc, "Stand", Date, "PYTHON", "")

except Exception, e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message

ERROR:

Line 18
ERROR 000539: : unexpected EOF while parsing (, line 1)
Fehler beim Ausführen von (CalculateField).


texts: Calculating the area covered by a polygon feature within a buffer area of a point in QGIS

texts: I have two datasets - one of points from bird nesting data and one of the extent of woodland. I have created 2km buffers around each nesting point and I want to find out what area of woodland there is within each point buffer. Can someone point me in the right direction as to what process I use here?


texts: database server connection to MSSQL Server with ArcCatalog not possible?

texts:
I have a MSSQL Server where my geodata is held. In the ArcGIS Docs, it says...

"To administer ArcSDE geodatabases for SQL Server Express, you can also connect to these geodatabases by first making a connection to a database server, then connecting to the geodatabase, as shown below."

But when I try it fails. Looking carefully again I notice that is applies to "SQL Server Express"....and not MSSQL Server. If I don´t have the express version of MSSQL does that mean I can´t create a database server connection? If I connect to an individual Database it works fine.

Thanks for any help,

Robert


texts: How to search for a specific field and then use it in Field Calculator

texts:
I have a piece of code where the goal is to search through the field names within a feature, select any field based on a specific character (in this case "z") and then use that selected field within a Field Calculator Expression. My experience with python is limited, see the code below;

for var in List_Input:

filename = os.path.splitext(os.path.basename(str(var)))[0]
arcpy.Select_analysis(var, env+"\\"+filename, "")
var_sel = env+"\\"+filename
arcpy.AddField_management(var_sel,"VR_Rate","Double")
VRField = arcpy.ListFields(var_sel,)

VRField_Calc = ""
for field in VRField:
    if "z" in field.name:
        VRField_Calc = field.name
        arcpy.CalculateField_management(var_sel,"VR_Rate", [VRField_Calc+" + 2"])

The bottom half of the code is where the issues arise. How to pass the "VRField_Calc" (field name) into the Calculate Field Expression?


texts: Exporting Map from Composer to PDF/PRINT, but all the work disapears

texts:
All my work in QGIS is missing when i tried to export via pdf, jpg or impress (including scales, legends, maps, layers, georeferenced photos). I'm not using any "OpenLayers Plugin".

I had already tried to export by "Maps Printer" Plugin. The legends and scale are there, but the map no.


texts: QGIS- Creating a python program to customize a login page for users of QGIS

texts:
Is there any method to create a customized login page for users to open QGIS desktop wherein the user activity can be tracked?
Once a user clicks on the QGIS desktop icon, the alert box shows login id and password and once the users gives the credentials . they can easily login to their particular project


texts: I am using QGIS 3.6 and I have a list of cities and zipcodes for my county, in my geopackage address layer I have a field called zip I tried to set up a value map with city names in the "description" and zipcodes in the "value" fields of the  widget that should allow me to select the name of the city the address is associated with and for the value I wanted it to place the zipcode in the field however when I select the description (or the city name) it places the name of the city in the zip value and not the zipcode IF I change the widget back to a text edit type it has the zipcode in the field. What is happening?


texts: Increasing performance without losing resolution when handling raster file in QGIS?

texts:
Orthomosaics generated by the processing of images obtained by UAV flights may often result in very heavy files that can take a long time to do simple tasks, just like to zoom on a specific region of the image.

I read about virtual raster and pyramid creation to deal with multiple raster files (on this topic -> Handling Many Raster Files in QGIS?) but I'm not sure if is this the best way to handle heavy raster files in QGIS?


texts: Rasterizing MULTIPLE vector layers in same "grid" so I can do a regression of sorts across pixel units?

texts:
For my research, I have collected a bunch of climatic data from my state in vector form. I also have my raster data, which is about 60 points around the state that have a binary assignment of yes/no. I'm trying to find out what climatic factors correlate with a point being either 'yes' or 'no'. My adviser suggested I rasterize my vector layers to 50m units and then start my analysis using that information.

My question is technical: is there a way to rasterize the entire map and collection of layers together, so that each point has an even chunk of each layer? Also, after rasterizing, can I export a CSV with the information generated for my analysis? (I know it would be huge)


texts: Translating gdal.polygonize from qgis 2 to 3

texts:
I was using this simple code to polygonize rasters using qgis 2.18:

general.runalg('gdalogr:polygonize', inRaster, 'DN', outShape)

I've tried converting it to use in qgis 3.2 like so:

general.run('gdal:polygonize', inRaster, 'DN', outShape)

But this error pops up and I can't seem to figure it out:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 16, in <module>
  File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\tools\general.py", line 96, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\core\Processing.py", line 132, in runAlgorithm
    feedback.reportError(msg)
AttributeError: 'str' object has no attribute 'reportError'


texts: Mercator: scale factor is changed along the meridians as a function of latitude?

texts:
"In Mercator projection the scale factor is changed along the
  meridians as a function of latitude"

How to calculate the Mercator scale factor as a function of the latitude in decimal degrees ?

I map any country of the world using it's bounding boxe's values.

1. In Equirectangular: The image width is linearly proportionate to the substraction of longitude values. The image height is linearly proportionate to the substraction of latitude values.

2. Moving to Mercator: The image width is linearly proportionate to the substraction of longitude values. The image height is linearly proportionate to the substraction of latitude values mutiplied by the scale factor at this average latitude.


texts: Re-positioning map grid coordinate labels in QGIS layout

texts:
I've created a map grid in QGIS (3.12) layout. Is there a way/expression to re-position the map grid coordinate labels so they are not overlying the map grid? It would be good to offset them to the side of the grids.

Is there also a way to add a drop shadow to these labels?

This is what the grid and labels currently look like:


texts: Get the whole coordinates of line string in OpenLayers API - ReactJs

texts:
I am working on a project where I am drawing a road with  OpenLayers map API as linestrings .
I could take the first and last coordinates of points of a the line.
This way:

drawend(e) {
    console.log('xxxxxxxxxxxxx, draw end', e);
    console.log("Co Ordinates :",e.feature.getGeometry().getCoordinates())
}

But I need to take data of all the line as an array.


texts: When iterating over a vector layer with pyqgis, how do I check whether a feature is selected?

texts:
While iterating over a vector layer using the following code (summerized from the example in the pyqgis cookbook), is there a way to check whether a feature is selected?

provider = vlayer.dataProvider()
feat = QgsFeature()
allAttrs = provider.attributeIndexes()
provider.select(allAttrs)
while provider.nextFeature(feat):
    geom = feat.geometry()
    attrs = feat.attributeMap()
    for (k,attr) in attrs.iteritems():
        print "%d: %s" % (k, attr.toString())

Alternatively, I could create a list of selected features using vlayer.selectedFeatures(), but I'm hoping there is a way to check each feature directly.


texts: Fill of sketch polygon blocks view of what I'm trying to digitize

texts:
I am trying to digitize polygons using QGIS 2.10.1-Pisa, but the red fill of the sketch polygon blocks me from seeing the imagery that I need in order to draw said polygon.

I have tried going into Settings > Options > Digitizing and changing the opaqity of the rubberband to a lower number (as described here: "add feature" tool QGIS 2.0 - polygon fill blocks view of what I want to digitize), but this does nothing. I believe it's a bug, but wanted to check to see if anyone else had figured this out.


texts: Connecting PostgreSQL/PostGIS and MapGuide Maestro

texts:
This is less a specific question, but hopefully more of an open question. There is little documentation regarding connecting the excellent OS web mapping software MapGuide Maestro (http://trac.osgeo.org/mapguide/wiki/maestro) with PostgreSQL/PostGIS (http://postgis.net/).

Can anyone give pointers on the connection of these two?


texts: Intersection between two ADF GraphicsLayer

texts:
I know how to intersect features from a ADF Geometry with a map ADF GraphicsLayer, but this requires me to have ONE LAYER on the map and then use a spatial filter on it.

If I have two GraphicsLayer which are NOT on the map, but they are the result of a query themselves, how can I intersect them and show the result?

I can't add one of them to the map because I can't interfere with it, this should be done in memory and not prompted to the user nor messing with the TableOfContents

Any idea on how to do that?


texts:
I am trying to use a conditional statement to generate a raster with binary values from a raster with probability values (floating point raster).  This is easily done within ArcMap by batching the process, but I need to do it >300 times with a different 'threshold' value each time.

Here is my code:

# Import arcpy module
    import arcpy

# Import environment and spatial analyst modules
    from arcpy import env
    from arcpy.sa import *
    arcpy.CheckOutExtension("Spatial")

# Set Workspace for Rasters
    env.workspace = "C:\May_2014\python_example.gdb"
    Rasters = arcpy.ListRasters("*", "All")

# Local Variables
    MyTable = "table_for_python"

    with arcpy.da.SearchCursor (MyTable, ["threshold"]) as mycursor:
        counter = 0
        for row in mycursor:
            raster = Rasters[counter]
            outCon = Con(raster >= mycursor, 1,0)
            outCon.save("C:\May_2014\python_example_2.gdb\Con_" + raster)
            counter = counter + 1

I keep getting the following error message after the line:

outCon.save('Con_' + raster)

AttributeError: 'int' object has no attribute 'save'.

Any suggestions would be most welcome.  Thanks in advance!


texts: Calculate number of alpha characters in field to a new integer field - ArcGIS

texts: In a feature class field (screenshot below), I'd like to know how I can get a count of the alpha characters in each record and calculate it to a new field.  I have a "FieldCount" field that calculates the total number of characters in the ENGPLANDWG field.  I'd like to do the same thing but just for alpha characters where "F63Y" would calculate to 2 (F & Y).  I'm using ArcGIS 10.2.2.


texts:
Not exactly sure how to word this, but hopefully I do it in a way that accurately conveys what it is that I'm trying to do...

I'm using Manifold GIS for this particular project, but I'm open to any methods that will achieve the desired result (less ArcGIS as I do not have any licensed ESRI software).

I started with a dataset that contained a series of points. Using a script posted by a Manifold user, I was able to convert these points into line segments (which may or may not be necessary to do what I need).

I have data like the following:

ID      US_NUMBER    Absolute_Distance_ft   DESCRIPTION    FEATURE_TYPE_ID   Line id
365694    5          0.059                  GWD            30                5
365695    0          0.059                  Weld Start     77                0
365696    0          1.814                  VLV 1          7                 0
365697    7          3.56                   GWD            30                7
365698    0          4.177                  Offtake        78                0
365699    10         4.783                 GWD            30                10
365700    0          5.774                  Weldolet       80                0
365701    20         6.709                 GWD            30                20
365702    0          7.602                                 5                 0
365703    0          9.209                  Bend           2                 0
365704    30         11.683                GWD            30                30
365705    40         20.705                GWD            30                40

Each line segment starts where the "FEATURE_TYPE_ID" is 30 and has a "US_NUMBER" of n. The next line starts at the next 30 and has a "US_NUMBER" n1, and so on...

As you can see in the table, "US_NUMBER" 5 is the start of one line segment. The next relevant line segment is "US_NUMBER" 7. The two segments between 5 and 7 are 0 and 0.

I want to go through the entire line and update each segments "Line ID" so that every segment that has a Line ID of 0 is set to the "US_NUMBER" where the previous "FEATURE_TYPE_ID" was 30...

So I would have something like...

ID      US_NUMBER    Absolute_Distance_ft   DESCRIPTION    FEATURE_TYPE_ID   Line id
365694    5          0.059                  GWD            30                5
365695    0          0.059                  Weld Start     77                5
365696    0          1.814                  VLV 1          7                 5
365697    7          3.56                   GWD            30                7
365698    0          4.177                  Offtake        78                7
365699    10         4.783                 GWD            30                 10
365700    0          5.774                  Weldolet       80                10
365701    20         6.709                 GWD            30                 20
365702    0          7.602                                 5                 20
365703    0          9.209                  Bend           2                 20
365704    30         11.683                GWD            30                 30
365705    40         20.705                GWD            30                 40

Does that make sense?


texts: Convert sfc-column into two columns with their corresponding x and y coordinates

texts:
I have a sfc-column with just 6 points. I'd like to make these point to markers in a leaflet-map in R. So I thought that it would be a nice idea to convert them into two columns: x and y. The sfc-column looks like the following. Probably an easy task, but I have no idea how to split up the values from each point into two columns.

Geometry set for 6 features
geometry type:  POINT
dimension:      XY
bbox:           xmin: 9.952508 ymin: 53.57191 xmax: 10.02894 ymax: 53.62756
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 5 geometries:
POINT (9.988762 53.5924)
POINT (9.975138 53.57954)
POINT (9.952508 53.57584)
POINT (10.02894 53.57191)
POINT (10.00046 53.59458)


texts:
I have a 3D map with latitude and longitude details. Is there a way I can take a 2D image (a png image) from 3D model in such a way that my 2D image also have the latitude and longitude as it is in 3d model ?

Can 2D image have the details as latitude and longitude ?

Which software should I use for doing this. Can anyone please point me to some online tutorial if available ?

Thanks,


texts: Projecting single pointGeometry object with ArcPy?

texts:
In an attempt to create a python add-in for ArcMap that will pass x and y coordinates to another application via a mouse click on the map, I have written the following code.  I need the output coordinates to be in WGS_1984 and therefore have to project a geometry object.  The problem I'm running into is that whenever the script executes the projectAs() method I get the following error.  Code at bottom

import arcpy
import pythonaddins
import webbrowser
import os, sys
from arcpy import env

arcpy.env.overwriteOutput = True

class ToolClass2(object):
    """Implementation for StreetView_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.

    def onMouseDownMap(self, x, y, button, shift):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        sr_df = df.spatialReference
        xValue = x
        yValue = y
        pointList = [[xValue, yValue]]
        point = arcpy.Point()
        for pt in pointList:
            point.X = pt[0]
            point.Y = pt[1]

            pointGeometry = arcpy.PointGeometry(point, sr_df)
            pointGeometry.projectAs("WGS 1984")

            centroid = str(pointGeometry.centroid)
            centList = centroid.split(" ")
            xValue2 = centList[0]
            yValue2 = centList[1]

I have tried adding the optional transformation for the projectAs() method as well without success.  What am I missing here?


texts: Can the XYZ tile naming convention be used for Geodetic (global-geodetic)?

texts:
I know that the XYZ convention is an alternative for the TMS convention for Mercator projection and in this configuration, it's also called the Slippy Maps naming convention.
Can the XYZ convention be used for the Geodetic projection? If so, is the conversion formula from TMS to XYZ the same for Geodetic and Mercator?


texts: Is there an existing object model diagram for ArcObjects? Something that can be printed out as a large poster.


texts: How I can extract a subset of years from a NetCDF file (e.g. to use in cell statistics)?

texts:
I want to use 5 years of data from a netCDF with 10 years of data (monthly slices). I saw this tutorial video that mentions this, but I can't figure out how to select more than one time slice.

http://video.arcgis.com/watch/2108/analyzing-netcdf-data-using-arcmap


texts: how to display a ee.Feature using folium?

texts:
I try to display a country delimitation in a jupyter notebook using folium.

I try the following code

#load country from LSIB
countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

#subset the sri-lanka
sri_lanka = ee.Feature(
    countries.filter(
        ee.Filter.eq('country_co', 'CE'))
    .first()
);

# Use folium to visualize the imagery.
mapid = image.getMapId(sri_lanka, {}, 'sri-lanka')

map = folium.Map(location=[80.773137, 7.873592])
folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; Google Earth Engine',
    overlay=True,
    name='border',
  ).add_to(map)

map.add_child(folium.LayerControl())
map

and I get the following error

  TypeError: getMapId() takes from 1 to 2 positional arguments but 4 were given

I don't really get the problem as in the documentation they describe it as:

  Get a Map ID for a given asset

  Returns the mapId call results, which may be passed to ee.data.getTileUrl or ui.Map.addLayer. Null if a callback is specified.

and this is exactly how I use it in gee.

Any idea ?


texts: Defining fields to keep and deleting all other fields in feature class using ArcPy?

texts:
I am trying to write a script to read a feature class, list all fields and keep only specic ones. I.e. drop all remaining columns.

I guess that pseudo code should be:

allFields = create list of all fields
fieldToKeep = define list of
fields to keep (manually)
fieldsToDrop = different items between
list 1 and list2 = to remove
execute
arcpy.DeleteField_management(inFC, fieldsToDrop) to remove
uncecessary fields

Should be pretty simple, but I cannot figure out steps and conditioning in correct order?

# Input point data
inFC= "inFC"

# CHeck if field exist,if yes, delete the field
allFields = arcpy.ListFields(inFC)

fieldsToKeep = ["OBJECTID", "Shape", "Country", "manualChange"]
fieldsToDrop = []  # define new empty list

for field in allFields:   # how to difference two lists?
    if not field in fieldsToKeep:
        fields.remove(field)

# Execute DeleteField
arcpy.DeleteField_management(inFC, fieldsToDrop)


texts: Implementing Add Feature action using PyQGIS

texts:
I'm trying to implement / reuse Add Feature functionality in my python plugin. there a class/function that handles geometry creation on canvas based on layer geometry definition? I'm checking API but I can't find a function that gives you functionality to digitize geometry and pass it as a feature to layer - or something similar.
I managed to get layer to editing mode now I can't find a way to implement this.
So basically --->  with PyQGIS API.


texts: I have lots of RADARSAT-1 file in GeoTiff format along with meta data where corner coordinates are available. Corner coordinates are for image portion(except the black portion). But, radar image has black shade on two sides, that make it difficult to georeferenced. How can I Geo reference all radar GeoTiff images. Can someone help me with code or some ideas?


texts: Is there a Feature-Feedback in QWC 2 (QGIS)?

texts: Is there a Feature-Feedback for QWC 2 (Qgis web client 2), where users can send some information for example wrong map information?


texts: How do I specify a spatial reference when using da.InsertCursor?

texts:
When using an InsertCursor, I can specify a spatial reference per
http://resources.arcgis.com/en/help/main/10.1/index.html#//018v0000002z000000.

I'm using da.InsertCursor because that's what everyone tells me to use and I've encountered no problems with it. How do I specify a spatial reference when using a da.InsertCursor?

Spatial reference is not even mentioned on the help page...
http://resources.arcgis.com/en/help/main/10.1/index.html#//018w0000000t000000


texts: I have two layers of Road Networks. How do I compare the difference in road direction between the two layers - both are polylines?


texts:
Below is code for calculating patch statistics on 34 rasters where each raster's name is written as the last column in the output. I'm trying to figure out what I'm missing to have the loop go through all rasters and write their results to one csv. Am missing an i somewhere or does it need to be a raster stack?

Update: when I turn list into a stack and the code below this is the error that I get: Error in if (m[ii] == 0) minp[ii] = 4 * n[ii] : missing value where TRUE/FALSE needed

library(raster) library(rgdal) library(SDMTools)
m <- "Z:/././"
allFiles <- list.files(m,pattern = "_ex0315.ovr", no.. = TRUE,recursive=TRUE)
allFiles<-stack(allFiles)

#PATCH STATS
install.packages("SDMTools")
library(SDMTools)

#Landscape patch analysis
#define a simple binary matrix

for (i in 1:length(allFiles)){
  A1 <- raster(allFiles[[i]])
  # plot(A1)
  A1[A1>0]<-1
  #do the connected component labelling
  ccl.mat = ConnCompLabel(A1)
  ps.data=list()
  #calculate the patch statistics
  ps.data[[i]] = PatchStat(ccl.mat)
  ps.data[[i]][1:nrow(ps.data[[i]]), ncol(ps.data[[i]])+1] <- allFiles[[i]]

}

Error in if (m[ii] == 0) minp[ii] = 4 * n[ii] :
  missing value where TRUE/FALSE needed


texts: CartoDB - update existing table from external data (i.e. Mysql database)

texts: I would like my CartoDB data to be updated, using external data (mysql database, I run and manage locally using php scripts), so that the CartDB table gets updated with the most recent data, whenever the map gets loaded. Is that even possible?


texts: Problems when trying to Undo edit operations in ArcGIS

texts:
I'm trying to programatically stop/abort/undo an edit operation inside a ArcSDE environment. My environment is currently Oracle, ArcSDE 9.3.1 and in  our custom code (I inherited a large ugly portion of a system) the developers do use StartEditing/StartEditOperation and StopEditing/StopEditOperation workflow, but they ONLY use IWorkspaceEdit interface and not Imultiversion interface.

Unfortunately, out database is NOT versioned and all edits done with ArcMap are done by uncheking the option of Undo/Redo - hence, all my undos fails (COM Exception).

I'm trying to capture a OnDeleteEvent and inside it, I check for a certain condition, and then I call AbortOperation, but the operation is not aborted. I need to stop the whole edit session and ask to not save to "rollback" my edits. This consistent with ESRI docs.

My guess here is, the delete already happened and commited to the database? If yes, what is an strategy to not allow this certain edit to happen? There are some other scenarios that I need to cancel the edit operation. I know that IEditEvents2 has a BeforeStopOperation event, but ESRI explicitely tells us not to call AbortOperation inside these events.

There are also other problems with this setup: I check for all edits on a stopEditOperation and generate some notifications with it. All edits come to me as empty. THis might have something to do with it?

Please, help?

Thanks guys!

Exemple code (onDeleteEvent):

    void events_OnDelete(IObject obj)
    {
        ITable table = null;
        IDataset dataset = obj.Table as IDataset;

        IFeatureWorkspace workspace = dataset.Workspace as IFeatureWorkspace;
        IWorkspaceEdit edit = workspace as IWorkspaceEdit;

        edit.EnableUndoRedo();

        IQueryFilter filter = new QueryFilterClass();

        try
        {
            foreach (string s in VariousTables())
            {
                table = workspace.OpenTable(s);

                if (s.Contains("vertex"))
                    filter.WhereClause = "clause A";
                else
                    filter.WhereClause = "clause B";

                if (table.RowCount(filter) > 0)
                {
                    edit.AbortEditOperation(); --> this still deletes the point
                    // edit.UndoEditOperation() --> this breaks, as the workspace does not support Undos
                    break;
                }
            }
        }
        catch (Exception ex)
        {
            MessageBox.Show("ErroR",ex);
        }
    }


texts: Print LABEL text in multiple lines in Cesium JS API

texts: I have a problem printing text in multiple lines using Label in Cesium JavaScript API. text is member of Label class which takes String as input but doesn't really support multiple line labels using '\n' character (typically as in Java string1 + '\n' + string2). Is there any way of printing labels in multiple lines?


texts: How to use UNION ALL in GeoServer CQL Filter?

texts:
Basically I want to combine 2 cql filters of the same layer and output the result on my geoserver layer. It is similar to this SQl query (my case using PostgreSQL):
select * from table1 where val = 'a'
union all  // union all means keeps duplicate
select * from table1 where val = 'b'

I tried using the cql_filter as below but didnt work
cql_filter="val"='a' union "val"='b'

I dont find the cql_filter union function (not union of geometries) in https://docs.geoserver.org/stable/en/user/filter/function_reference.html as well. What is the proper way to use UNION ALL likewise the SQL query in GeoServer cql filter?


texts: The resampled image has many NaN values and does not cover the complete area of the input raster

texts: when I resample an image of 20m to 240m resolution, the resampled image does not cover the whole area of the input image and does the resampling for some pixels in which I get many NaN values. How can I solve this problem?


texts: Using field calculator to change domain type

texts:
I have a table with a bunch of domains, and a table with the information needed to populate the domains.  I joined the two, but the table of information is populated as text.

I need to be able to change the domains to what the information table has.

I have checked and the domains and the information table have the same number of selections.

Is there a way in the field calculator using python to select and change the coded value?

so if I have 3 choices :

shapefile domains:

1 A

2 B

3 C

information table:

A

B

C

Is there a script that would change A to 1 B to 2 C to 3 and have the domains change?


texts: LAStools in QGIS error at activation: "wrong value for parameter 'OTB applications folder'..."

texts:
I have a problem with LAStools in QGIS 2.14. I've followed the LAStools instructions for downloading LAStools for QGIS. `

I replaced
'lidar' folder in C:\OSGeo4W64\apps\qgis\python\plugins\processing\algs\lidar with the folder downloaded from http://rapidlasso.com.

I downloaded the 'LAStools' folder and put it in the following location: 'C:\LAStools'. I activate LAStools in to Processing, and type 'C:\LAStools' for the 'LAStools folder location', but I get the following message:

  "Wrong value for parameter "OTB applications folder":
  Specified path does not exist:
  C:/OSGEO4~1/apps\orfeotoolbox\applications"

LAStools appears in the processing toolbox. I can start to run a lasview. The data appears in 'just a little LAS and LAZ viewer', but the program hangs and then crashes.

Is this happening because of the error message? I have looked through the StackExchange for other examples of this specific problem, with no success. Any suggestions on what am I doing wrong?


texts: Obtaining elevation data for points using QGIS

texts:
I am new to QGIS.

I have points uploaded to QGIS that were collected from a GPS, and I need to get elevations for each point (I don't trust the elevation estimations of my GPS).

What is the best way to do this?

Additionally, I realize I'll likely need an elevation raster, and I'm wondering if the Google terrain hybrid or Google satellite hybrid will work or if there's a better option?


texts: QGIS atlas, add picture with path from attribute

texts:
I love working with the atlas generator and I came up with a new case to use it for. I want to add a picture to my map from URL.
I am allowed to add a picture if I type in the full url path. eg. http://my.url.service.no/pics/img_001.jpg
But I want to add the picture based on the filename from my table.
See example in my print screen.

Using the same expression creating a lable it give me the string for the complete url eg. http://my.url.service/[%"filename"%] evaluates to  http://my.url.service/img_001.jpg

Is it possible to do what I try here, anyone with experience or tips?


texts: What's a term for Sea Level Distance which is usable by novice users?

texts:
I've been using the term Sea Level Distance to describe the total length of a line string where I ignore the elevation and just calculate the distance from point to point along the ellipsoid and sum the total.

Sea Level Distance has always seemed to be an artificial term for this. Slope Distance seems to naturally describe the distance when elevation is included. Is there a term for Sea Level Distance that is more natural and usable for novice users.


texts:
When I type gdal_merge.py into Terminal, the output is -bash: gdal_merge.py: command not found.

I installed 3.4.0 Madeira [from this link] for macOS 10.13 (https://www.qgis.org/en/site/forusers/download), which came with a gdal installer.

I cannot even find the gdal_merge.py file in /Library/Frameworks/GDAL.framework/Versions/Current/Programs/ or anywhere else in the gdal framework directory. I looked to see if I could install the .py file separately but can't find any installation options on the GDAL website.


texts:
I need to extend several lines to its intersection with another feature class (polygons).

In order to know the distance from the road to the building facade (in several points of the road) I created perpendicular lines to the street segments each 5 meters and I want to extend the line, until the building facade.

I have tried the 'Extend tool' in QGIS but it only let me extend them to the buildings if they are in the same layer (if they are lines and then merging the building lines with the perpendicular lines). However, I do not consider it as the best method because if all lines (building lines and perpendicular lines) are in the same layer perpendicular lines could just intersect themselves (and not the building lines).

Is there any way in ArcGIS Desktop I can make those lines extend until they intersect with the buildings as a separated feature class?

Is there any way I can make it in QGIS?

If any of you know a script or something I could also try using Python.


texts: Multiple Leaflet TileLayers with extended maxzoom

texts:
I've got 2 layers(for indoor and outdoor) in my web application. If i zoom beyond 19, the indoor layer should be displayed and the outdoor layer tiles will be scaled up, because their maxzoom is 19. This works fine with this js-code:

map = L.map('full-map');
tileLayer = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 23,
    maxNativeZoom: 19
}).addTo(map);

var indoorLayer = L.tileLayer("http://" + SERVER_IP + ":8002/tiles/" + 0 + "/{z}/{x}/{y}.png", {
    minZoom: 20,
    maxZoom: 23
}).addTo(map);

Except this is how it looks when i add the indoor layer to the map:

What makes me wonder is how some outdoor tiles are still displayed.

It works perfectly if i don't add the indoor layer to the map:


texts: How to create multible labels non-ovelaped in CartoDB?

texts:
I'm trying to make a map pointing the cities I've visited for work, in each city I've visited more than one place. Those places are my labels, bur they are overlapped. If I choose label overlap = false, only one label shows, but I need all of them to appear.

Can anyone help me?


texts: Silence user warnings for GRASS-PYTHON program

texts:
I've created a program that runs externally from GRASS v6.4.1, using python 2.6. It imports an ESRI shapefile, and exports a seperate KML for every polygon in this vector file.

Polygons with islands produce a warning message using "v.out.ogr" saying that  features without a category were skipped.  Looking into this further, these "skipped" features (islands) are processed later.  So the program does what it was designed to do, but I would like it so that this warning message is not printed in the terminal.

VERBOSITY setting were set to "0" using:

os.environ['GRASS_VERBOSE']='0' #runs modules silently

But this doesn't take care of warning messages.  I found that there is a gisenv variable called DEBUG, can be set to "0" as well so that all warning messages are silenced too.  Problem is, I can't seem to get the proper syntax for my GRASS-python script.  I've tried:

grass.gisenv() ['set']='DEBUG=0'

Using this line of code atleast allows the program to run, but warning messages are still printed.  Other options I've tried terminate the program.

Has anyone successfully set the DEBUG level to "0" in a GRASS-Python script?


texts: Filtering GPS tracks to determine "Moving Time"?

texts:
I am working with GPS tracks of outdoor activities, think cycling, running, hiking. These activities often have "rest periods" where people take a break, wait at traffic lights, or just stopped their GPS for whatever reason. I end up with the raw GPS track and I am wondering how to best calculate "moving time". All gis.stackexchange.com I have looked at are very vague, and I am wondering if there is literature/methods for doing this.

I can of course go the easy route and simply look at all trackpoints and subtract all those time periods where the speed between to trackpoints is under some threshold. I would need to tweak this threshold for the activity type (I am assuming), but this seems not very elegant. I have various other ideas, but it seems to me that this problem must be common, and maybe there's a well-known solution?

Also, if anybody has any info how services such as Strava are doing it (as simple as I described above, or something more sophisticated?), it would be great to get some info about that as well.


texts: Merging separate polygons to one polygon without using (buffer, dissolve, and neg buffer)?

texts:
I want to merge separate polygons to one polygon. But I want to do it without using (buffer, dissolve tool ets). I want to visualize the sepated polygons in to one polygon.

I can work in ArcGIS and Qgis.