texts: Serialize geometry to GeoJSON with ArcObjects

texts:
In the ArcObjects developer SDK, I found that the IExternalSerializerGdb2 interface of the JSONSerializerGdb class supports writing to GeoJSON.

  ArcGIS Developer Help (ESRI.ArcGIS.Geodatabase)

  IExternalSerializerGdb2.WriteRecordSetEx2 Method

  Writes feature
  recordset for both JSON and GeoJSON format and returns number of
  records serialized.

I tried this:

var jsonWriter = new JSONWriterClass();
jsonWriter.WriteToString();

var jsonSerializer = new JSONSerializerGdbClass();
jsonSerializer.InitSerializer(jsonWriter, null);

var recordset = new RecordSetClass();
recordset.SetSourceTable((ITable)featureclass, null);

var sr = ((IGeoDataset)featureclass).SpatialReference;

int count;
((IExternalSerializerGdb2)jsonSerializer).WriteRecordSetEx2(null, recordset, null, sr, out count);

But it only seems to write the (ArcGIS) JSON format.

How do I generate GeoJSON?


texts: Adding features from ArcGIS Online to QGIS?

texts:
How can I add features which provided by ArcGIS Online to QGIS map?

I'd like to add feature like "USA Mean Temperature" with QGIS (2.18.3).
These useful features are easy to access when I use ArcGIS Online,
by doing "Add" -> "Search for Layers" -> Select "ArcGIS Online".

I guess "Layer" -> "Add Layer" -> "Add ArcGIS FeatureServer Layer" would be answer, but I'm not sure what URL I should use.

Could you tell me how I should do?


texts: Creating 3D terrain from point data in PostgreSQL?

texts:
I am quite new to arcgis and I would like to create a 3d terrain. Now I have a postgresSQL table containing a list of points with their geometry data and the Z value for each point.

What I would like is guidance on how to proceed in order to create this terrain in 3D in arcgis.

I might be taking the wrong approach so please can someone guide me on how I can create this 3D terrain from points data that I have taken with my gps?


texts:
I got an old vector dataset with polygons covering a continent. The data was first published on paper at scale 1:5 000 000 and was later digitalized. I don't have the original data and no information about the vectorization or any metadata. I guess that the distance between the vertices rather than the accuracy limits the resolution.

Vertices are saved with high resolution (e.g. "nnn.nnnnnnnnn","-nn.nnnnnnnnn"). The dataset has few points that can be georeferenced nor any nodes that are defined as coordinates (e.g. at even degrees or UTM coordinates). When I compare some coastline sections, the error is up to +/- 20km.

I'd like to find a formula to estimate the maximum error based on the distribution of the vertices. I have access to any GIS application but would prefer a robust statistical reference.

How can I calculate the maximum error of the dataset, assuming that all vertices are correct? Or phrased differently: What method can I use to find the largest resolution of the dataset?

I tried to rasterize the dataset at different cell sizes and then oversample it to a small cell size to detect the smallest possible rasterizing without loss of resolution, but that is rather time-consuming and not very mathematical approach.


texts: Get intersection in polygons with holes or multi parts

texts:
Like in a title. Let's say I have normal polygon that overlaps other. And in my case, when the other polygon has holes or is splitted (like on attached image), normal methods "intersection" or "intersects" don't work.

One of my try was something like this, but it failed.

"some_feature" is that normal polygon.

# if polygon has holes, create new one without them
if len(feature.geometry().asPolygon()) > 1:
    feat = QgsFeature()
    geom = feature.geometry().asPolygon()
    feat.setGeometry(QgsGeometry.fromPolygon([geom[0]]))
    print(some_feature.geometry().intersection(feat.geometry()).area())

Below are some images, that show cases I'm talking about.

On the left is case with holes and on the right is one object as two polygons.

Does anyone have idea, how to achieve such thing?


texts:
I am having a hard time putting labels on both ends of a line
Here is SLD that I have so far. It puts a label on one side only but I need it on both ends of a line.

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
        xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
        xmlns="http://www.opengis.net/sld"
        xmlns:ogc="http://www.opengis.net/ogc"
        xmlns:xlink="http://www.w3.org/1999/xlink"
        xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
    <NamedLayer>
        <Name>line_name_ends_red</Name>
        <UserStyle>
            <Title>line_name_ends_red</Title>
            <FeatureTypeStyle>
                <Rule>
                    <Name>line_name_ends_red</Name>
                    <LineSymbolizer>
                        <Stroke>
                            <CssParameter name="stroke">#ff0000</CssParameter>
                            <CssParameter name="stroke-width">2</CssParameter>
                        </Stroke>
                    </LineSymbolizer>
                    <TextSymbolizer>
                        <Geometry>
                            <ogc:Function name="endPoint">
                                <ogc:PropertyName>the_geom</ogc:PropertyName>
                            </ogc:Function>
                        </Geometry>
                      <LabelPlacement>
                       <PointPlacement>
                         <Rotation>
                          <ogc:Function name="endAngle">
                                    <ogc:PropertyName>the_geom</ogc:PropertyName>
                                </ogc:Function>
                         </Rotation>
                       </PointPlacement>
                      </LabelPlacement>
                      <Label>
                         <ogc:PropertyName>LINE_NAME</ogc:PropertyName>
                       </Label>
                    </TextSymbolizer>
                </Rule>
            </FeatureTypeStyle>
        </UserStyle>
    </NamedLayer>
</StyledLayerDescriptor>

But I need it to look like this


texts:
Based on Creating sectors for telecom towers?, I'm facing one problem after creating the site sectors.

I can sucessufuly create all sectors with shape tool. After that I save my project and everything is ok.

But if I close QGIS and reopen the project, my shape layer that should have all sectors, is now empty without features on attribute table. (image below)

Can anyone explain me why this is happening?


texts:
I am trying to create a base map suitable for overlaying data that is in SRS EPSG:4326 using CartoDB following the projections tutorial.  I have other projections working, but not one for this.

So using:

  SELECT cartodb_id, ST_Transform(the_geom, 4326) AS the_geom_webmercator
  FROM world_borders_hd

It does not seem to return any geometries plotted in the tiles.

From a comment on the CartoDB academy issue I posted, I presume this is something to do with a scale that needs to be applied in the transformation from degrees to meters.

Can anyone advise me what I need to do in the SQL to get a zoom 0 tile that would look like this?


texts: Bulk changing annotation locations in ArcMap?

texts:
I'm working with a large number of labels that I've converted to annotation (using ArcMap 10.1). All of them need to fit within the parcel borders that they correspond to.

Is there a way to do it all at once or am I stuck going one by one to fit them in the boundaries?


texts: Will core functions in QGIS exploit multi-threading?

texts:
After reading Enabling multithreaded processing in QGIS?, I was wondering if QGIS 2.6 will include this when running various processes. I checked online and I could only find that multi-threaded rendering was introduced to QGIS 2.4 (I assume it will again be available in 2.6). Now maps can be redrawn very quickly especially when dealing with large datasets.

I read that it's fairly difficult to include multi-threading to core functionality and that users are advised to manipulate Python code in order to achieve this, such as in Parallelising GIS operations in PyQGIS?.

I also checked the feature request but it has been closed for 9 months now by Tim Sutton with the last comment being:

  "I'm closing this out - Martin Dobias has an implementation in a branch which will be merged post QGIS 2.2"

Will some functions in QGIS 2.6 exploit multi-threading (or will it again be focused on rendering) and if not, will that change at QGIS 3.0?


texts: Programmatically adding shapefile to data frame in ArcMap?

texts:
I have a toolbox with python script which creates a shapefile from JSON txt file, but it does not show the generated shapefile in arcmap, i have to add it from the output location.

What code shall I use to display this shapefile once it is created?


texts:
I'm using the semi-automatic classification (SCP) plugin in QGIS on a Mac OS Catalina 10.15.5. I've tried it in multiple versions and I'm having issues with saving files. Whenever I click the Run button on a Mac, it automatically creates temp files rather than letting me save the outputs somewhere.
For instance, when trying to clip multiple rasters using a PC, I click the Run button and a window pops up allowing me to browse to a folder to save the files. On a Mac, a window opens allowing me to only click Open rather than save and the output is a temp file with no discernible file name to help me figure out which of the multiple clipped files it is referencing and doesn't retain the default 'clip' prefix.
Any ideas? I can't find anything online that references how to fix this.


texts: Finding if latitude and longitude in shapefile region

texts:
I referenced this very helpful question, but still cannot crack the issue.

I downloaded Zillow's neighborhood shapefile for Connecticut and have a latitude and longitude: 41.344602, -73.071443 which exists in CT, according to Google Maps.

As far as I can see, the shapefile is using the same format, yet when I use the following snippet:

import shapefile
import shapely.geometry as geometry
shapefile_location = 'C:/Downloads/zillow_ct_neighborhood_boundaries/ZillowNeighborhoods-CT.shp'
shp = shapefile.Reader(shapefile_location)
pt = 41.344602, -73.071443
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(0,len(all_shapes)):

    boundary = all_shapes[i] # get a boundary polygon
    print( all_records[i][2],boundary.points)
    name_ =  all_records[i][2]
    if 'New Haven' in name_:
        ok_regions.append(boundary)
    #     print(name_,geometry.Point(pt).within(geometry.shape(boundary)))
        if geometry.Point(pt).within(geometry.shape(boundary)): # make a point and see if it's in the polygon
            name_ = all_records[i][2] # get the second field of the corresponding record
            print("The point is in", name  )

I also tried to reverse the latitude, longitude combo, but no luck. My code does not find my latitude, longitude combo in any of the neighborhoods.

Does anyone have some guidance on what I may be doing wrong?


texts: Problem downloading and using NOAA CCAP data

texts:
I am trying to download CCAP land cover data and would have expected .img or .tiff files, but the data downloads as a disk image file. I cannot figure out how to access/work with this file type. any help?

I want to use the data in R to create time series of land cover change in different basins.

Here is a link to download the data: https://coast.noaa.gov/ccapftp/#/
I am attempting to download the Virginia Regional Dates datasets


texts: Finding if field has empty records using ArcPy

texts:
I am trying to iterate all fields in a feature class/shape file to check for empty values. If there are empty records in any field, I would like to print that field.
I know Python a little bit and here is my attempt.
import arcpy
from arcpy import env

env.workspace = r"C:\SOFTWARE\Python\Jun2019"
fc ="test.shp"
##
fieldlist = arcpy.ListFields(fc)
##Iterate feilds
for field in fieldlist:
           # Iterate over each feature; if it has a non-null value
           with arcpy.da.SearchCursor(fc, field) as cursor:
               for row in cursor:
                   if row[0]h=

I am stuck here.
How do I continue?

I want print fields TYPE,CAT,SUBTYPE etc. which are having empty values


texts: How to Save line/polygon location data as csv for later import

texts:
I have some Line layers that I would like to save as as a .csv file, edit the attibute columns and the re-import to be saved as a shapefile.

How can I add geometry or location attributes to be saved as .csv, and how do I correctly manage them when importing? Would also be interested to know if this can be done for Polygons.


texts:
I would like to integrate the time slider into a JSONP geoserver. Here is the current code I have http://jsfiddle.net/qcwcrwge/26/ and I would like to use the time slider with a layer called feature.properties.time.

I have tried to add the time slide libraries and the lines

var testlayer = L.geoJson(json),
        sliderControl = L.control.sliderControl({
            position: "topright",
            layer: feature.properties.time
        });

//For a Range-Slider use the range property:
sliderControl = L.control.sliderControl({
    position: "topright",
    layer: feature.properties.time,
    range: true
});

//Make sure to add the slider to the map ;-)
myMap.addControl(sliderControl);
//And initialize the slider
sliderControl.startSlider();

Just after the line

}).addTo(coolPlaces);

But it does not work.

Do you know what the problem can be?


texts: How to generate heatmaps with QGIS1.9 plugin automatically?

texts:
I created a few heatmaps with the new plugin in QGIS 1.9, but now I would like to generate a lot of them automatically. How should I do?

I created a few of them by hand, but I should first import the data, select the CRS, then select the plugin, select the imported data and the parameters for the heatmap (and pick again the CRS). It takes quite long as soon as you want more than a few heatmaps.

I have 300'000 data in one day. I separate them in groups of 15 minutes, and I would like to visualize the heatmap for all groups of 15 minutes (~64), in order to observe the evolution in time.

I'm interested in either:

a way to make it in QGIS automatically, or
a way to do it in another software, that would be more adapted for this purpose.


texts:
I have two line datasets. Both the datasets are falling in the same region. Now few of them are matching spatially with each other when I viewed them from top. I wanted to identify all the lines which are matching with each other (both partially and fully). None of them are having common attributes.

How can I do this?


texts: Optimal method to convert raster from 32bit Float to 8bit Byte to reduce file size in QGIS 3

texts:
I have a raster file (.tif) that contains integer values from 1 to 31. I have realized that the values are stored as Float32 with NoData Value=-3.39999999999999996e+38 (see raster information below). I think that the optimal method to store the data is as unsigned 8 bit integer (QGIS Byte) as there are no decimal values and the maximum value is below 255. I have used the Translate (convert format) tool available from GDAL to change the data type from Float32 to Byte (following these instructions). However, I have run into two issues:

The resulting file is BIGGER than my original file (from 31MB to 268MB). I noticed that the original file had COMPRESSION=LZW and tried to use various compression formats under  Profile (i.e. default, no compression, Low compression, High compress, JPEG compression) but I was not successful in reducing file size while maintaining data quality (e.g. using JPEGwith 100 resulted in some values being +/-1 of the values they should be and some random additional pixels)
I had difficulty setting the NoData value. It seems that 0 worked in my case but I may have files with 0-values that are meaningful and I wouldn't know how to convert the multidecimal NoData value from the Float to an integer value that is not confused with actual zero values in the dataset.

Metadata of original file

Size is 21768, 12920
Pixel Size = (250.000000000000000,-250.000000000000000)
            Image Structure Metadata:
              COMPRESSION=LZW
                  NoData Value=-3.39999999999999996e+38

UPDATE:
According to this related post, the raster size is computed based on the following formula:

  The size of a raster is just the product of bit-depth/8, bands, rows,
  and columns plus header metadata (statistics, etc.)

32/8*   1*  21768*  12920*  =1124970240 Bytes/(1024*1024) ~1072MB
8/8*    1*  21768*  12920*  =281242560  Bytes/(1024*1024) ~268MB

That gives me the size Windows indicates for my 8-bit raster without compression.

These are the options in the dialog box:


texts: Image equalization after tiling, using gdal_translate -scale

texts:
I want to generate a set of samples from a Sentinel-2 image by tiling it, using mainly gdal_translate. The tiling is ok, but its histogram seems to vary tile to tile. So, what I thought was to apply the same min/max from original image to each tile, but didn't workout.
After to search and use examples found here and here, the tiles still present variations. The code used:
name, file_extension = os.path.splitext(image)
image_path = os.path.join(image_folder, image)

ds = gdal.Open(image_path)

stats = [ds.GetRasterBand(i + 1).GetStatistics(True, True) for i in range(ds.RasterCount)]
        vmin, vmax, vmean, vstd = zip(*stats)

rows = ds.RasterXSize
cols = ds.RasterYSize
tiles_cols = cols / width
tiles_rows = rows / height
logging.info(">>>> Tiling image {}. {} x {} pixels. Estimated {} tiles of {} x {}..."
             .format(image, rows, cols, round(tiles_cols * tiles_rows), width, height))

gdal.UseExceptions()
for i in range(0, rows, width):
    for j in range(0, cols, height):
        try:
           output = os.path.join(output_folder, name + "_" + str(i) + "_" + str(j) + file_extension)
           gdal.Translate(output, ds, format='GTIFF', srcWin=[i, j, width, height],
                          outputType=gdal.GDT_UInt16, scaleParams=[[list(zip(*[vmin, vmax]))]], options=['-epo'])
        except RuntimeError as error:
           logging.info(">>>>>> Partially outside the image. {} Not saved!".format(error))

when I open the tiles individually, it seems to have different histograms, even the stretching parameters are unique (calculated from original image, see stats) to every generated tile.

Am I missing something or should I see a uniform tiled-image?


texts: Can I upload an LPA file (Long Profile) to QGIS

texts: I have a .lpa file that I'd like to upload to QGIS - can this be done and if so, how?


texts: Display two different maps with different layers in one layout in PyQGIS? Problem with .pdf export

texts:
I would like to export a .pdf with 2 maps (with different layers).
However, some layers do not appear in the map items!

I'd like:

In the first map: the shapefile over the raster 1, and these two over the basemap (A300 + Area1 + basemap)
In the second map: the shapefile over the raster 2, and these two over the basemap (A300 + Area2 + basemap)

Everything appears fine in the layout view (see image left) but in the .pdf only the basemap appears (see image right)! What is wrong? here is my code

import os
import processing
from qgis.core import *
from qgis.gui import *
import qgis.utils
import shutil
from qgis.PyQt import QtGui

project = QgsProject.instance()
project.setCrs(QgsCoordinateReferenceSystem(21781))

iface.mainWindow().blockSignals(True)

# add basemap (raster)
basemap = "C:\\MyData\\LK25.vrt"
basemap_r = iface.addRasterLayer(basemap , "my_basemap")
basemap_r.setCrs(QgsCoordinateReferenceSystem(2056))

# add vector layer
shplayer = iface.addVectorLayer("C:\\MyData\\A300.shp", "A300", "ogr")
shplayer.setCrs(QgsCoordinateReferenceSystem(2056))
project.addMapLayer(shplayer)

# add raster layer 1
rlayer1 = iface.addRasterLayer("C:\\MyData\\Area_1.asc", "Area1")
rlayer1.setCrs(QgsCoordinateReferenceSystem(2056))
project.addMapLayer(rlayer1)

# add raster layer 2
rlayer2 = iface.addRasterLayer("C:\\MyData\\Area_2.asc", "Area2")
rlayer2.setCrs(QgsCoordinateReferenceSystem(2056))
project.addMapLayer(rlayer2)

#-- Create custom caneva
manager = project.layoutManager()
layout_name = 'my_custom_layout'
layout_list = manager.printLayouts()

layout = QgsPrintLayout(project)
layout.setName(layout_name)
layout.initializeDefaults()
manager.addLayout(layout)

#--- create map item 1 (shapefile, raster 1, basemap)

map1 = QgsLayoutItemMap(layout) # create Layout item
map1.setRect(QRectF(10, 10, 210, 90)) # set the map item rectangle
map1.setPos(10,10) # map item upper left corner position
map1.setFrameEnabled(True) # to have a map black border

map1.setLayers([shplayer, rlayer1, basemap_r]) # set layers to be mapped
map1_extent = QgsRectangle(rlayer1.extent())
map1_extent.scale(1.0)
map1.setExtent(map1_extent)

layout.addLayoutItem(map1)
map1.attemptMove(QgsLayoutPoint(10,10))
map1.attemptResize(QgsLayoutSize(210,90))

#--- create map item 2 (shapefile, raster 2, basemap)

map2 = QgsLayoutItemMap(layout)
map2.setRect(QRectF(10, 110, 210, 90))
map2.setPos(10,110)
map2.setFrameEnabled(True)

map2.setLayers([shplayer, rlayer1, basemap_r])
map2.setExtent(map1_extent)

layout.addLayoutItem(map2)
map2.attemptMove(QgsLayoutPoint(10,110, QgsUnitTypes.LayoutMillimeters))
map2.attemptResize(QgsLayoutSize(210,90,QgsUnitTypes.LayoutMillimeters))

layout.refresh()
project.write()

# export map in .pdf
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C:\\qgs2.pdf", QgsLayoutExporter.PdfExportSettings())


texts: OpenStreetMap and Country Borders in QGIS?

texts: I am new to QGIS and I am using OpenStreetMap and OSM place search to create a map for latin america. However the peruvian terrritory includes the maritime area. How can I avoid/modify this. All other countries are fine.


texts: Creating Aeronautical maps covering entire Colombo Flight Information Region (FIR)?

texts:
ICAO need Aeronautical maps in WGS84 Coordinate System but in ArcGIS length can't be given or measure as it's a GCS.

In Sri Lanka we have SLD99(EPSG5233), Kandawala/Sri Lanka Grid(EPSG5234) and UTM44(EPSG32644) PCSs.

What PCS or PCSs should be used to create Aeronautical maps (Aerodrome, Approach and En-route) for the whole Colombo FIR?


texts:
I'm new to Webmapping I already create a map with ReactLeaflet and ReactLeaflet Draw, but I can't imagine how can I get the geoJSON from the polygon and just send to a database or an API
    export default class Draw extends Component {

  _onEdited = (e) => {
    let numEdited = 0;
    e.layers.eachLayer( (layer) => {
      numEdited += 1;
    });
    console.log(`_onEdited: edited ${numEdited} layers`, e);

    this._onChange();
  }

  onCreated = (e) => {
    let type = e.layerType;
    let layer = e.layer;

    if (type === 'marker') {
      // Do marker specific actions
      console.log("_onCreated: marker created", e);
    }
    else {
      console.log("_onCreated: something else created:", type, e, );
      alert('Foi criada uma layer!!');
      console.log('criada com sucesso');
    }
    //Here I need to Stringify the Layer into GeoJSON for use it later.

  }

  _onDeleted = (e) => {

    let numDeleted = 0;
    e.layers.eachLayer( (layer) => {
      numDeleted += 1;
    });
    console.log(`onDeleted: removed ${numDeleted} layers`, e);

    this._onChange();
  }

  render() {
    return (
    <>
      <Map
      className="leaflet-div"
      center={[37.8189, -122.4786]} zoom={13} zoomControl={false}>
        <TileLayer
          attribution='&copy; OpenStreetMap contributors'
          url="http://{s}.tile.osm.org/{z}/{x}/{y}.png"
          id= ''
        />
        <FeatureGroup>
            <EditControl
                position="topright"

                onEdited={this._onEdited}
              onCreated={this.onCreated}
              onDeleted={this._onDeleted}
                edit={{ remove: false }}
                draw={{
                    marker: false,
                    circle: false,
                    rectangle: false,
                    polyline: false
                }}
            />
        </FeatureGroup>
      </Map>
//When the user click on the button I need to send the info to an API with AXIOS.
      <button>
     Send Polygon
    </button>
      </>
    )
    ;
  }
}


texts: How to improve performance using spatial indexing with STIntersects()

texts:
I am working on Spatial database. I need to check intersection of state geometry with geometry inside database. I am presently using the query as:

SELECT
    a.ID,
    a.OBJECTID,
    a.MSID,
    a.FACNAME,
    Convert(varchar(max),a.ShapeGeo),
    b.State_Name
FROM
    XFAcilityBackUp26June a with(index(testIndexSpatial)),
    us_stateTest b
WHERE
    a.shapegeo.STIntersects(b.the_geom.MakeValid()) = 1 and b.State_Name = 'Texas'

It is functioning properly but it takes more time(more than 10 minutes) to execute. I have tried indexing also but no luck with it.

Is there any alternative way to improve the performance?


texts: Making network /tower coverage map in Mapbox?

texts:
I am new to mapping and currently I am trying to make a map showing coverage of a network provider here in my country. The map I am trying to achieve should look similar to that at http://www.imore.com/best-us-iphone-carrier-2014.

I already have the latitude and longitude coordinates of the tower locations in CSV format and have so far managed to plot the markers in TileMill (ideally I want to show coverage as opposed to markers). I have also downloaded Mapbox studio and OSGEo4W under the assumption that I might have to convert the data into another format.

I can't seem to find any comprehensive tutorials for this, all I see for Mapbox is OpenStreetMaps, terrain, markers, shapefiles etc.

I have also attached an image of what I am trying to achieve


texts: Calculate and plot the feature importance of the input dataset of a random forest classification within Google Earth Engine (GEE)

texts:
I run a random forest classification for agricultural land use and other land cover classes (12 classes). My dataset for 2019 consists of a Sentinel 1 and Sentinel 2 monthly time series, statistical phenometrics (e.g. min/ max/ median/ stdDev) and ancillary data as texture metrics and topographic derivates. I included a lot of different indices calculated on the basis of the sentinel 2 optic bands. This resulted in a single image with 294 bands as a big input data cube for the random forest algorithm. In order to decrease computational time I would like to calculate the feature importance of the different bands in that image for each one class against all classification and then only use the most important bands. So far I didn't find a lot in different blogs and platforms.
I found an approach on MyGeoBlog that is unfortunately showing an error that I couldn't resolve. In theory this approach should work with the one class against all other classification approach. Find my adjusted script with the following link (note that its a different dataset with only 216 bands without the sentinel 1 data):
https://code.earthengine.google.com/?scriptPath=users%2Fmarcellitomohr93%2FArmenia_LULCC%3ALULC_Classification2019%2FGoris2019_ClassificationIterations%2FFeatureImportance_GrasslandMask
The original script on MyGeoBlog you may access via this link (however it is lacking input data):
https://mygeoblog.com/2020/02/12/random-forest-trees-and-variable-importance/
In my case GEE shows the usual point feature size exceeding error but also the error: "Error generating chart: Property 'Grass_Mask' of feature '0' is missing." Since for this script you can only calculate probability with 2 classes I transformed my 12 classes of a LULC_Class property to a binary property column that is called grass mask (with only values of 1 for grassland and 2 for all other classes).
Im happy to hear about your thoughts and experiences with calculating the feature importance in the GEE environment. I could also do it in R but then I have to download my very large dataset.
With best regards,
Marcel


texts: QGIS Raster reprojection error (gdalwarp: no such file or directory)

texts:
I am working with QGIS 3.4.4 Madeira (OSGeo4W)
I have a Digital Elevation Model/Digital Terrain Model of a watershed with a EPSG:3006 - SWEREF99 TM - CRS (that is actually unsuitable for the DEM cause the DEM-Area is out of the boundaries of that CRS. That file was given to me by someone else, so I don't know why they used this crs)
I wanted to reproject it to WGS84 UTM33N. I used the QGIS Raster->Projection->Reproject it returns the error ERROR 4: false: No such file or directory

that is the command (i replaced the file path wit ...)

{ 'DATA_TYPE' : 0, 'EXTRA' : False, 'INPUT' : 'C:/...file.tif', 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'C:/.../file.tif', 'RESAMPLING' : 0, 'SOURCE_CRS' : None, 'TARGET_CRS' : QgsCoordinateReferenceSystem('EPSG:32633'), 'TARGET_EXTENT' : None, 'TARGET_EXTENT_CRS' : None, 'TARGET_RESOLUTION' : None }

GDAL command:
gdalwarp -t_srs EPSG:32633 -r near -of GTiff false C:/Users/Tobias/Documents/Studienprojekt/DGM/Lockwitz_DGM.tif C:/.../file.tif
GDAL command output:
ERROR 4: false: No such file or directory

The file is there! I tried with the file loaded to QGIS; directly put in the file-path; save as temporary-file; save permanent
Nothing works always this Error, I also tried with another DEM (ASTER-GLOBAL-DEM, and SRTM-DEM) but I always get this error. So is there something I am doing wrong or misunderstand regarding that function? Or has anybody experienced similar issues with that? (Just started using QGIS so excuse if I overlooked something)


texts:
I have a simple code that is pulling together polygon shapefiles from a Windows folder and merging them together into a single shapefile. In the same code, I want to add a dissolve command. My fields are "Name" and "Status". If there are two polygons with the same Name (i.e. 2 different rows with the same Name), I want to dissolve by Name to create unique rows of data.

I am using ArcGIS 10.1. Please see my code below:

import arcpy, os, time

arcpy.env.workspace = r'C:\\Temp\\Planning\\Planning_In\\Validated'
fcs = arcpy.ListFeatureClasses(feature_type='Polygon')

fcOut = r'C:\\Temp\\Planning\\Planning_Out\\ProjectMerge.shp'
if fcOut in fcs:fcs.remove(fcOut)

if len(fcs) > 0:
    if arcpy.Exists(fcOut):
        fcArchived = os.path.join(r'C:\\Temp\\Planning\\Planning_Out\\Archive', 'Archived_{}.shp'.format(time.strftime('%Y%m%d')))
        os.makedirs(os.path.dirname(fcArchived))
        arcpy.Copy_management(fcOut, fcArchived)
        arcpy.Delete_management(fcOut)

print fcs

arcpy.Merge_management(fcs, fcOut)

# Set local variables
inFeatures = r'C:\\Temp\\Planning\\Planning_Out\\ProjectMerge.shp'
tempLayer = r'C:\\Temp\\Planning\\Planning_Out\\tempLayer.shp'
expression = arcpy.AddFieldDelimiters(inFeatures, "Name") + " <> ''"
outFeatureClass = r'C:\\Temp\\Planning\\Planning_Out\\ProjectDissolve.shp'
dissolveFields = ["Name"]

# Execute MakeFeatureLayer and SelectLayerByAttribute.  This is only to exclude
#  features that are not desired in the output.
arcpy.MakeFeatureLayer_management(inFeatures, tempLayer)
arcpy.SelectLayerByAttribute_management(tempLayer, "NEW_SELECTION", expression)

# Execute Dissolve using LANDUSE and TAXCODE as Dissolve Fields
arcpy.Dissolve_management(tempLayer, outFeatureClass, dissolveFields, "",
                          "SINGLE_PART", "DISSOLVE_LINES")

I end up with a shapefile, but it is not dissolved. Also, my "Status" field that is in the existing, singular shapefiles is not carried over.

The statuses of all projects will be the same. For example, if there are two polygons with Name= 1415-006, they will both have the same status before the dissolve takes place. If I use the following code and also dissolve by "Status", I have the "Status" field in the dissolve shapefile but still, the items are not dissolved.

# Set local variables
inFeatures = r'C:\\Temp\\Planning\\Planning_Out\\ProjectMerge.shp'
tempLayer = r'C:\\Temp\\Planning\\Planning_Out\\tempLayer.shp'
expression = arcpy.AddFieldDelimiters(inFeatures, "Name") + " <> ''"
outFeatureClass = r'C:\\Temp\\Planning\\Planning_Out\\ProjectDissolve.shp'
dissolveFields = ["Name", "Status"]

# Execute MakeFeatureLayer and SelectLayerByAttribute.  This is only to exclude
#  features that are not desired in the output.
arcpy.MakeFeatureLayer_management(inFeatures, tempLayer)
arcpy.SelectLayerByAttribute_management(tempLayer, "NEW_SELECTION", expression)

# Execute Dissolve using LANDUSE and TAXCODE as Dissolve Fields
arcpy.Dissolve_management(tempLayer, outFeatureClass, dissolveFields, "",
                          "SINGLE_PART", "DISSOLVE_LINES")


texts: Can you get the individual segment length from OpenLayers2 measure class?

texts:
OpenLayers has a OpenLayers.Control.Measure class to implement a measure tool.

There is a measure and measure-partial event, both which are passed a variable containing the total distance. Is it possible to get the last segment length?

CLARIFY: When I say last segment length... you see how you can draw 2 or 3 lines in one go and it will give you the total distance drawn across all lines? That's great. What I'm after is the length of the last line the user drew only.


texts: weighting and combining data sets in Raster Calculator

texts:
I`m having problems using Raster Calculator to create a final cost raster. I have to add 7 raster data sets together and multiply each of them with the correct weight. The resulting raster is completely empty with unrealistic values of High:3,40282e+038 and Low:-3,40282e+038.

All the layers have the same cell size, same extent and no NoData. Any ideas what might be going wrong there?


texts: Add feature id as a feature property in Google Earth Engine

texts:
In GEE, user-generated feature collections are automatically assigned an id number. How do I create a property for each feature collection that stores the id ("system:index") as a property, such that when I export the feature collection as a shapefile, the ID remains?
I'm trying to execute this with the variable tiles_new before exporting the feature collection. Code is here.


texts: Randomly coloring provinces in country map?

texts:
I have created a map of a country in QGIS.

How do I color the province individually ?
Is there a quick way to randomly coloring all the province of a country.

I am using the "ne_10m_admin_1_states_provinces" shapefile from Natural Earth.


texts: Drawing ellipse filling dataframe with ArcPy?

texts: I need to modify one of my template scripts to draw a circle from the centre of the (inset) data frame to the edge of the data frame. This will then be used to clip the data frame to shape. I'm using ArcGIS 10.0.


texts: ArcGIS field calculator: VBScript get data from previous row

texts:
When using the field calculator, how do I get data from the previous row in VBScript? In Python I'd do something like this:

for index, row in enumerate(rows,start=1):
    row[index-1] += 5

I'm using ArcGIS 9.3.1.

Let's assume I have a column [Distance] and all it's rows are 0. How do I increment every row with 5 starting at row 1. So the first row remains 0, second row'll be 5, third 10 and so on.


texts:
Good day all,

I am attempting to understand if (and how) I might persuade QGIS to perform a union operation on the drawing of the marker symbols (not the geometry itself) based on a series of points.  It's easy to set up SVG markers, but let's say that what I see after doing so is this:

What I'd like to see, however, is this:

I know this is easily accomplished if the green circles were actual polygon geometries, but what I'm wanting united (unioned?) are the SVG marker symbols as drawn on the screen/printed - not the point data.

Is this possible, if so, how?


texts: Selecting objects from a feature layer by attribute using OR?

texts:
I need to select all objects from a feature layer where a certain field called ESTATE meets 'rule a' OR 'rule b'. I have managed to do so with the following piece of code:

arcpy.MakeFeatureLayer_management(fc_Optimity_All, 'fl_Optimity_All')
query = "ESTATE LIKE '{}'".format('Node Live')
arcpy.SelectLayerByAttribute_management('fl_Optimity_All', 'NEW_SELECTION', query)
query = "ESTATE LIKE '{}'".format('Node Planned')
arcpy.SelectLayerByAttribute_management('fl_Optimity_All', 'ADD_TO_SELECTION', query)

What I have done above is basically an OR query (using the ADD_TO_SELECTION option on the second select instruction). Is there any way to do the same thing but only with one "arcpy.SelectLayerByAttribute_management" instruction?


texts: How to change text color in QGIS 3 composer label when render as HTML

texts:
I try to set a label with black background, a first line of text in green (exactly rgb(125,185,0) ) and the other line in white.

In QGIS 2.18 setting the background black, the font color white and using HTML tag on the first line <Font color = rgb(125,185,0)> worked.

On QGIS 3.6 the first line color change but to a kind of blue (see picture) and trying other RGB combination produce color that doesn't follow any RBG code, it seem that the color tag has an effect but that QGIS doesn't use RGB to render it (in QGIS color selector the rgb(125,185,0) produce the expected green).

Is that a QGIS 3 bug or are my setting wrong?


texts:
Ever since google maps came out, the standard for web maps has been in their “flawed” web Mercator projection. I am trying to figure out why this projection is not good from a geospatial analysis perspective. Here is what I am trying to understand:

Say I create a buffer around a point that is 1 mile wide. This calculation was done in ArcGIS using WGS84 as my coordinate system. So, I will have a lat/long for all of the points around my single point within a mile. If I then want to show this on a google map, it must be projected into web mercator, which uses meters as a coordinate system. If I trace this reprojected 1 mile buffer, how closely will it resemble my buffer that used a projection using lat/long instead of meters? Will it run over all of the same houses/trees/etc if I am looking at the imagery basemap?

Or, if I have “google coordinate system/projection” on Arcmap and perform my buffer process on this, will it be less accurate than if I had done it on WGS84? Does reprojecting between the two result in a noticeable difference? Can analysis be done in a meter coordinate system?

Sorry if my questions are a bit vague. I just want to figure out what types of problems I may incur performing analysis and having the results on a google map. Thanks a lot!


texts:
I have used QGIS to batch import and visualise tiled elevation data. This is illuminated greyscale, with identical parameters used each time, so I would have expected the result to be seamless and the tile boundaries to be essentially invisible.

However, the result is a 'patchwork quilt' effect. The shading is not consistent from one tile to the next. I do not understand why not. This is the code used to apply the hillshading:

  def ShadeRaster(raster, pth):

    import processing

    parameters = {'INPUT': "",
      'BAND': 1,
      'COMPUTE_EDGES': False,
      'ZEVENBERGEN': False,
      'Z_FACTOR': 1.0,
      'SCALE': 1.0,
      'AZIMUTH': 315,
      'COMBINED': False,
      'ALTITUDE': 45,
      'MULTIDIRECTIONAL': False,
      'OUTPUT': ""}

    parameters['INPUT'] =raster.name()
    parameters['OUTPUT'] = pth + '/' + raster.name() + '.tif'

    processing.runAndLoadResults('gdal:hillshade',parameters,feedback=feedback)
    return QgsProject.instance().mapLayersByName('Hillshade')[0]

As you can see, the elevation and azimuth of the light source are always the same. But this is how it looks:

Why?


texts: Determine the number of fields of a certain field type there are in an input feature class

texts:
So I am wondering how to get the number of field types "String" by creating a custom function? Can anybody point me in the right direction? I have a shapefile named streets and there are 73 field names/ field types and I keep getting that number and not the corresponding number for only type "String".

import arcpy
arcpy.env.workspace = "C:\\Users\\Ogi\\Desktop\\GIS Programming\\Data\\PrgGIS\\Lab10"

def countstringfields(table):
    fields = arcpy.ListFields("streets.shp", {"String"})
    namelist = []
    for field in fields:
        namelist.append(field.type)
    return len(namelist)

fieldtype = countstringfields("C:\\Users\\Ogi\\Desktop\\GIS Programming\\Data\\PrgGIS\\Lab10\\streets.shp")
print fieldtype


texts: Interrupting contour lines beneath elevation labels (instead of using label buffers) using QGIS?

texts: Is there a way to interrupt a contour line beneath an elevation label using QGIS?


texts: Writing only final raster with intermediate bands going to disk using ArcPy?

texts:
I have a raster with 4 bands.

I would like to create a new raster (its individual bands called new_band) from this raster (its individual bands called old_band) in with :

new_band1 = old_band1 * 2
new_band2= old_band2 * 4
new_band3= old_band3 * 6
new_band4 = old_band4 * 8

My code is:

import sys, os, math, time
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

raster = r"D:\IMAGERY.tif"
bands = [['band1', 2], ['band2',4],['band3',6],  ['band4', 8]]
outfolder =r"D:\Radiance"

for band in bands:
      num = band[0][4:]
      thisBand = raster + '\\band_' + num
      Radiance = arcpy.sa.Float(Raster(thisBand)/float(band[1]))
      outname = 'B' + num + '.tif'
      Radiance.save(os.path.join(outfolder,outname))
#and then combine different bands
  arcpy.CompositeBands_management("D:\Radiance\B1.tif;D:\Radiance\B2.tif;D:\Radiance\B3.tif;D:\Radiance\B4.tif", "D:\CANGIO\Anh_SOPT5\Radiance\compbands.tif")

In this code, I created both images for individual bands and a new raster.

As I just would like to create a new raster file, not rasters of individual bands, how can I do this?


texts: In an attempt to summarise landclasses, the columns of my pivot table output are a bit mixed up with the field info regarding the three input fields that I've selected. As shown in snapshot 3 the description field have the info of the grid_slo field, the grid_slo field have the info of the grid_alt field, and grid-alt showing zero. The original descriptions that should show under the descriptio field has been lost. Numerous attempts to fix this problem has been unsuccessful. Can anyone help please?


texts:
I am attempting a GIS tutorial on land change analysis between two rasters. I have set the extent and raster analysis to what I have been given and used the expression below. However, it is not working. Could someone provide some clues as to what I have be doing wrong? I have followed all instructions... I can't see what the error no is.

Con("%lc1991%" == "%lc2001%",0,1)

000539 is the error I get and I had a look and it said cannot concatenate a string and integer. This is the expression I was told to use. Does anyone know how to construct a better expression?


texts: How can I get the latitude and longitude bounds of a tile set from MapTiler using leaflet.vectorgrid?

texts:
I have a leaflet map that loads a vector tile set from Maptiler. I would like the map to fit the bounds of the tile set's content when the user selects the overlay check box to display this tile set. For example, I'm working with the zip codes in Illinois. Each zip code is a polygon.

How can I get the minimum and maximum latitude and longitude of this tile set? The tiles are saved in a protobuf (protocol buffer). If this is impossible, is it at least possible to get the x and y coordinates of the centroid of a single tile? Or at least the tile's z, x and y indices? Then I can use those values to calculate latitude and longitude.

I am using the leaflet.vectorGrid library.

  var ILZipTileOptions = {
    interactive: true,
    vectorTileLayerStyles: {
      ILZipCodes: properties => {
        //maybe I can add a line here that gets the coordinates?
        return getColor(
          [properties.ZCTA5CE10]
        );
      }
    }
  };

  var ILZips = L.vectorGrid
    .protobuf(IL_ZIP_CODES_PBF, ILZipTileOptions)
    .on("click", function(e) {
      L.popup()
        .setContent(
          "Zip Code: " +
            e.layer.properties.<zip code variable>
        )
        .setLatLng(e.latlng)
        .openOn(myMap);
    });

After poking around some more, I see it is possible to get the x, y and z indices of the tiles as they load.

  var ILZips = L.vectorGrid
    .protobuf(IL_ZIP_CODES_PBF, ILZipTileOptions)
    .on("load", function(e) {
      console.log("e.target._tiles: " + JSON.stringify(e.target._tiles));
    });

Here is an example printed out in the console...

e.target._tiles: {"3:6:4":{"el":{"_leaflet_pos":{"x":153,"y":200}},"coords":{"x":3,"y":6,"z":4},"current":true,"retain":true,"loaded":1579817242542,"active":true},"4:6:4":.....

But I cannot yet figure out how to programatically grab the x and y values from this object.

Update:
I figured out how to programatically grab the x and y values.
The way I'm accessing the variable that has colons in the variable name is ugly. Maybe there is a better way.

  var ILZips = L.vectorGrid
    .protobuf(IL_ZIP_CODES_PBF, ILZipTileOptions)
    .on("load", function(e) {
      console.log(
        "tile_x_coordinate: " +
          e.target._tiles[Object.keys(e.target._tiles)[0]].coords.x
      );
      console.log(
        "tile_y_coordinate: " +
          e.target._tiles[Object.keys(e.target._tiles)[0]].coords.y
      );
      console.log(
        "tile_z_coordinate: " +
          e.target._tiles[Object.keys(e.target._tiles)[0]].coords.z
      );
      console.log("e.target._tileZoom: " + e.target._tileZoom);
    });

The next step I took was to convert these tile indices to latitude and longitude degrees like this...

  function tileIndicesToDegrees(tile_x, tile_y, zoomLevel) {
    var n = Math.pow(2, zoomLevel);
    var lon_deg = (tile_x / n) * 360.0 - 180.0;
    var lat_rad = Math.atan(Math.sinh(Math.PI * (1 - (2 * tile_y) / n)));
    var lat_deg = (lat_rad * 180.0) / Math.PI;
    return [lon_deg, lat_deg];
  }

At the end of all this work I realized that the MapTiler tile set covers the entire world and not strictly the area of interest (i.e. the area covered by the geojson I uploaded to MapTiler). Therefore, this task of using the tile set indices to find the boundaries of the original geojson features is pointless.

How can I get the boundaries of the geojson features I used to create this tile set? Or is this impossible?


texts: How to set color of single features with Python

texts:
I have a layer with cities. I want to set a different color for a specific country. I can iterate through all cites (features) and change the color of ALL features but not only of single ones. How can I achieve that?

from PyQt4.QtCore import *
from PyQt4.QtGui import *

#Get layer
layer=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "cities_min":
        layer = lyr
        break

#Iterate features. How to set color of single features
for f in layer.getFeatures():
    print f['name'],f['country_code']
    if f['country_code']=='GB':
        #Set Color
        break

#Apply color to ALL features
symbols = layer.rendererV2().symbols()
symbol = symbols[0]
symbol.setColor(QColor('yellow'))
layer.triggerRepaint()


texts: Export CSV from CartoDB with lat and lon?

texts: If you create a point dataset in CartoDB then export it to CSV, you get a "the_geom" column in the WKB format. How do you get individual lat and lon columns?


texts: How to sum up values within grid polygon in QGIS?

texts:
I have a grid full of hexagons that I am trying to bin points. However, instead of using the Vector > Analysis Tools > Points in Polygon which counts the number of point inside each hexagon, I'm trying to sum up the values of other attributes for each hexagon. Better yet, it would be great if I could do this as well as count the points inside each hexagon.

For example, my datatset is formatted like this:

ID  Lat  Long  Score
1  43.1 -431.5 .75
2  44.3 -431.6 .89
3  40.8 -429.7 .12
4  100.1  80.9 .36
5  99.8   76.1 .65
6  91.2   75.9 .24

What I'm hoping will be something like this:

Hexagon Sum   N
  1     1.76  3
  2     1.25  3


texts: I have seen many examples on leaflet of data being plotted onto a google maps style map, is it possible to plot leaflet data onto a map with only coastline? What I'm looking for is just a white background with a black line tracing the coast, nothing else.


texts: How to compute time-space maps automatically?

texts:
I am looking for a way to produce a time-space map in QGis or ArcGIS. See an example here.

As input data, I have a city in the centre of Germany with travel times to some other cities around it. I am looking for a way to deform the "normal" geographic map of Germany to show the travel time instead of distance.

Any idea on how to do that with QGIS, ArcGIS or other GIS ?


texts:
I am new to both python and QGIS and I am trying to load a shapefile into python and I keep getting an error
when I try to read it.

sf=shapefile.Reader('wetland6.shp')

ShapefileException                        Traceback (most recent call last)
<ipython-input-48-dc4b8aeb80b5> in <module>()
----> 1 sf=shapefile.Reader('wetland6.shp')

C:...\shapefile.pyc in __init__(self, *args, **kwargs)
    232         if len(args) > 0:
    233             if is_string(args[0]):
--> 234                 self.load(args[0])
    235                 return
    236         if "shp" in kwargs.keys():

C:...shapefile.pyc in load(self, shapefile)
    265                 self.shp = open("%s.shp" % shapeName, "rb")
    266             except IOError:
--> 267                 raise ShapefileException("Unable to open %s.shp" % shapeName)
    268             try:
    269                 self.shx = open("%s.shx" % shapeName, "rb")

ShapefileException: Unable to open wetland6.sh


texts: Why does QGIS fail to include authentication manager overlays when using standalone exportToPDF?

texts:
I have an existing layout and I'm trying to export it to PDF from a standalone script after the data is updated automatically. If I open QGIS and export all is good but when I call it from Python I get the basemap and point overlay only. My polygon layer is missing.
I had to unlock/lock the layers then save the layout to get the points to display in the exported PDF but no joy for the polygons.
Both the point and polygon layer come from an Oracle database layer but only the polygon layer uses the authentication manager for the credentials (i.e. authcfg).
Code example:
#!/usr/bin/env python3
import os, time
from qgis.PyQt.QtCore import QSettings
from qgis.PyQt.QtSql import QSqlDatabase, QSqlQuery
from qgis.core import (QgsProject, QgsLayoutExporter, QgsApplication, QgsAuthManager, QgsFeedback)

QgsApplication.setPrefixPath("/usr", True)

gui_flag = False
app = QgsApplication([], gui_flag)

app.initQgis()

project_path = os.getcwd() + '/project.qgs'

project_instance = QgsProject.instance()
project_instance.setFileName(project_path)
project_instance.read()

authMgr = QgsApplication.authManager()
if authMgr.authenticationDatabasePath():
    print('INFO: Auth already initilised => we are inside a QGIS app.')
    msg = 'Inside QGIS'
else:
    print('WARNING: outside qgis, e.g. in a testing environment => setup env var before')
    # db init
    os.environ['QGIS_AUTH_DB_DIR_PATH'] = "profile/qgis-auth.db"
    msg = 'Master password could not be set'
    assert authMgr.setMasterPassword("password", True), msg
    authMgr.init( "profile/qgis-auth.db" )

if not QSqlDatabase.isDriverAvailable('QOCISPATIAL'):
    print("WARNING: Oracle not available")
else:
    db = QSqlDatabase.addDatabase("QOCISPATIAL")
    db.setDatabaseName("oracleDatabase:1521/dbName")
    db.setUserName("username")
    db.setPassword("password")

    if db.isValid():
        if not db.open():
            print("WARNING: Oracle not available")
        else:
            print("INFO: Oracle database available")

manager = QgsProject.instance().layoutManager()
timestr = time.strftime("%Y%m%d-%H%M%S")

feedback = QgsFeedback()

layout = manager.layoutByName("layoutName") # name of the layout
# or layout = manager.layouts()[0] # first layout

exportSettings = QgsLayoutExporter.PdfExportSettings()
exportSettings.flags = layout.renderContext().flags()
exportSettings.dpi = 150
exportSettings.forceVectorOutput = True
exportSettings.writeGeoPdf = True
exportSettings.includeGeoPdfFeatures = False

exporter = QgsLayoutExporter(layout)
exporter.exportToPdf(project_instance.absolutePath() + "/layoutName-" + timestr + ".pdf",
                     exportSettings)

app.exitQgis()

You can see that I have already put in a number of, largely unnecessary, checks to ensure that it is connecting to the Oracle database and those seem to be working. In addition it now draws the point layer which is also coming from the Oracle database.
My working hypothesis is that the authentication credentials are failing but I don't get an error so it is difficult to debug that. Since changing from QgsLayoutExporter.PdfExportSettings() to defining them manually I now get errors about the Dimensions of the inserted PDF page differing from the output PDF page but I'm not sure they are linked as the PDF is still created successfully; just without the polygon layer.
Update
I have now tried the script again using the username and password typed into the project and I can confirm the layer is always missing if I use an authcfg and works fine if the username and password are in the project file.
My question is therefore changed to; how can you use the authentication manager for a layer and still be able to export to PDF using a standalone python script?


texts: How to get a copy of the raw photo from the Mapillary app

texts:
In Mapillary's documentation it says the photos are stored in the "app folder" before upload, and are deleted afterwards.  I have looked around the filesystem but can't find them.  I don't necessarily want to have to upload a photo to get a copy of it; and after upload, the raw photo isn't available anyway - it's already been overlaid with a "Mapillary" watermark.

details: I am using the Android client.


texts: A script to populate a numeric field in Arcgis based on some content of a string field?

texts: I have field in an attribute table populated with expressions like "6S2Olch1E", "10S1B", "9B3S2Olch" and so on. I need to make a numeric filed populated with only numbers associeted with letter S from those experessions (f.e. 6 if expression is "6S2Olch1E" or 10 if it's "10S1B"). Is there a way to do this using a script in Field Calculator? I'm not a programmer.


texts:
I'm trying to create a points in the middle of polylines. The data is now in the form of shapefile.
Is there an easy way to do this?


texts: ArcMap: in a 10.0 Basic license and trying to edit FC coming from a 10.0 Editor edited GDB

texts:
I have a feature class in a file geodatabase which is not allowing me to edit it, but is also not giving an error when I try.  I right click on the layer in ArcMap and choose to start an edit session, Arcmap starts one and if i look in the Editor menu it does show it is in an editor session, but the layer I am trying to edit shows it is NOT open for editing (ie. "Start Editing" is not greyed out).  There is no error message saying i am not licensed to edit that data.

I have had this before when working on data sets previously worked on in a higher Editor/Advanced license but I don't know what part of the puzzle keeps the dataset from being edited by lower license levels.  There is no linked Annotation or relationship class type functionality to the dataset, just a domain on it and that's it, as far as i can tell.

...I pulled the FC into a blank MXD and tried to edit and this time DID get an error dialog that i was not licensed to edit the data, so I think that IS what's happening but for some reason editing it in the other MXD does not bring up the error dialog (other MXD likely opened in the higher license previously). So what do I need to do in the higher license to make the lower licenses able to edit the dataset?

The end user is on 10.0 and not upgrading anytime soon, so I have to make this backwards-compatible as well as license-agnostic.


texts: Drawing two points in two datums with same coordinate value

texts:
Today I was trying to use ArcMap to find out what is the distance between two locations with same coordinate values in NAD27 and NAD83, but the points I drew appeared to be at the exact same location. There must be something wrong with my understanding in datum or my knowledge in ArcMap, but I can't figure it out. Below is what I did:

Open a blank mxd
Create a Feature Class A in a File Geodatabase in NAD27 coordinate system
Create a point in Feature Class A at (40, -75), save edit
Open a new mxd
Create a Feature Class B in a File Geodatabase in NAD83 coordinate system
Create a point in Feature Class B at (40, -75), save edit
Load Feature Class B in the first mxd

The result was that the points in two feature class appeared at the exact same location. I believe there is something to do with projection on the fly, but I thought if the two points represented two different real world locations, they should appear in different locations on a map no matter what coordinate system the map uses. What am I missing here?


texts: How to build a Spatialite road network from OSM data

texts:
I'm looking for algorithm which could help me with building my own road network (for routing) from OSM data.

I chose Spatialite DB for this, but here I have only one problem: Spatialite provides interface for build road network, but prepared road graph is saved in BLOBS so I can't get any information about it (of course, I can write queries to DB for routing facilities) but I want to draw all links on my map, so I need to get their coordinates.

Actually, I'm not sure if Spatialite could provide such API or smith like that.

My question is how can I build my own road graph (network) from raw OSM data (I've already parsed XML file for highways=*) ?


texts: How to batch label features using PyQGIS?

texts:
I am writing a python script to label features (lines). However I can't figure out:

how to get the script to turn on the buffer
how to get the script to apply the changes. When I currently run it, it changes everything in the properties (except the buffer) but I still have to hit ok before QGIS applies it to the map.

.

layer = iface.activeLayer()
layer.setCustomProperty("labeling/fieldName", "Length" )
layer.setCustomProperty("labeling/placement", QgsPalLayerSettings.Horizontal)
layer.setCustomProperty("labeling/predefinedPointPosition", QgsPalLayerSettings.BottomRight)
layer.setCustomProperty("labeling/fontSize","10" )
layer.setCustomProperty("labeling/buffer",QgsPalLayerSettings.BufferDraw)
layer.setCustomProperty("labeling/enabled","true" )
layer.triggerRepaint()

I have had a look at one similar question but I cannot seem to figure it out:
How to label line features using python code?

Once I get this to work I want to transform it into a batch-script by adding some lines like:

for layer in QgsMapLayerRegistry.instance().mapLayers().values():


texts: Is pre-1990 historical block-level census data available online?

texts:
I'm looking for historical block-level data for Chicago.

According to this  Census bureau document, historical block-level census data exists for major cities going back to 1960, and even back to 1940 (though only for housing statistics.)

Only 1990 - 2010 is available on the census website.

Do digital copies of historical census block boundaries and stats exist? Are they available anywhere?

If not, where are the analogue versions of this data stored? And what is the most granular level of historical census data that is available digitally?


texts:
I'm looking to get QGIS Server working on a Mac with El Capitan OS.  When I try to access this link:

http://www.qgis.org/wiki/QGIS_Server_Tutorial

it gets forwarded to this link:

https://issues.qgis.org/projects/qgis/wiki/QGIS_Server_Tutorial

So I created an account for the QGIS issues tracker, but I still can't access the file. I get error "403. You are not authorized to access this page"

Can someone help me access the QGIS Server Tutorial or explain what I need to do to get the server working on El Capitan?


texts: How to separate zones into individual maps

texts:
I am using ArcMap 10.3 with a Basic license.

I have a single map split into four zones.

I would like to split this single map into four individual maps / data frames focusing on a single zone.

Is this possible? I know that I can just zoom to an extent and print that, but I would like to completely separate the content of Zone A from anything else.

In the above example, I wouldn't want to see the data on the edge of Zone D in my Zone A map.

I'd like to avoid duplicating and renaming my data four different times in order to just highlight and delete the excess data.


texts: Seeking software to visualize and export area selected in OSM format?

texts: I have a very large OSM file and I would like to know if there is software that can open the OSM file and then when I select a certain area to be able to export it in OSM file again. I found Maperitive and Marble but both of them can't export area in OSM format.


texts: Is it possible to link an Odoo database (CRM) into GGIS? I have an online Odoo Studio account with lots of data including long/lat data etc. Is it possible to link the Odoo database directly with Qgis or is that limited to PostGis etc?


texts:
Is there a way to get the extent (bbox) of a vector layer in QGIS?

I see that I can update the extent, but I am looking for the actually coordinates of the extent.


texts: Joining table with thousands of transactions to one point (location) in ArcMap?

texts:
I have a large transaction data set in the millions and around a 100 positional points. I need to join the transaction data set to these points. ArcMap won't let me join them without errors or if it does allow it null values are applied to the join.

They have a common field, a meter number for the machines. I have tried using a simple join with the attribute tables. My end goal is to be able to upload this data to an enterprise and integrate it into an operations dashboard. This dashboard would be able to tell the user the number of transactions with the machine within a time period. The transaction data set also has the date attached to it.

I am using ArcGIS Desktop 10.0.

I found a problem. The one web service map I am using has the meters as strings so in the transactions attribute table I added a new field and made it a text. Using the field calculator I added the meter number I was then finally able to perform the join, but only 153 values are in this join meaning it didn't join everything I wanted.

How can I get all of the values to join?


texts:
I am trying to export a line feature class from an oracle spatial database to a file geodatabase.  I can export polygons with no problem, but cannot find the appropriate command '-nlt' for Line Feature.

    ogr2ogr  -update -overwrite -nln  CONVERT_TO_LINEAR  -nln sf_line -a_srs EPSG:3089 -progress -f FileGDB g:/...." "


texts: Checking value of field based on conditional if statement using ArcPy cursors

texts:
I'm using ArcGIS desktop 10.6 and executing an arcpy script.
I'm creating a script that checks an input shapefiles field values against a set of rules. To do so I've created 2 lists with unique values of 2 fields within my input shapefile.
These are the UniqueDispTypeValues list and RelUnqID list.
I also have created a list of allowable values to check against.
This is the DispTypeValues list.
I created an if else statement to check a basic rule (If only RTF value exists the Rel_Unq_ID must be 0), however I want to check the current value of the field after the if statement(kind of validating a field value, or passing a rule in my case) instead of just the print statement.
Here is a screenshot of the attribute table:

Here is the code:
# Importing necessary modules
import arcpy, os

arcpy.env.overwriteOutput = True

# Setting path to mxd and Data Frame
mxd = arcpy.mapping.MapDocument(r"CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

# Assign the OneSTOP shapefile to a variable
shp = r"P:\JOBS\2020\20-1095\Outgoing_Information\OneSTOP\Disposition.shp"

#Create the unique list of allowable Disposition Types
def DTV(CHK_TABLE , Disp_Type):
    with arcpy.da.SearchCursor(CHK_TABLE, [Disp_Type]) as cursor:
        return sorted({row[0] for row in cursor})

DispTypeValues = DTV(r'Q:\GIS\OneSTOP\OneSTOP2.gdb\Checklist' , 'Disp_Type')

print (DispTypeValues)

#Create the unique list of Disposition Types in the input shapefile
def UDTV(shp , Disp_Type):
    with arcpy.da.SearchCursor(shp, [Disp_Type]) as cursor:
        return sorted({row[0] for row in cursor})

UniqueDispTypeValues = UDTV(shp , 'Disp_Type')

print (UniqueDispTypeValues)

#Create the unique list of Rel_Unq_ID in the input shapefile
def RUID(shp , Rel_Unq_ID):
    with arcpy.da.SearchCursor(shp, [Rel_Unq_ID]) as cursor:
        return sorted({row[0] for row in cursor})

RelUnqID = UDTV(shp , 'Rel_Unq_ID')

print (RelUnqID)

length = len(UniqueDispTypeValues)

print (length)

for value in UniqueDispTypeValues:
    if value == 'RTF' and length == 1:
        print "Rel_Unq_ID must be 0"
    else:
        print "Rel_Unq_ID must be non zero"


texts: Is the undo functionality possible in QGIS while creating a line feature?

texts:
I have noticed that when creating a new line feature in QGIS, control+z does not seem to be possible. When you insert a vertice of this line in a slightly wrong place, you cannot undo your actions (remove vertice). Control+z is only possible after finishing the feature, and in that case the whole feature is removed, and not the last verticed added.

Anything that solves this problem?


texts: Delete Pseudo Nodes (or unncessary nodes) using PostGIS or QGIS

texts:
I have so big big data (like about 60 Gb.) and it seems to have much unnecessary pseudo nodes.

I would like to remove them using PostGIS or QGIS without simplifying or changing its original shape
(Simplification might helpful though it changes the original shape and this is not what I want... also, ArMap DataReviewer is not going to be helpful or GT Tools such as ArcMap Addin)

What function names should I look into?
I have polygons and polylines that need to be worked done.

Finding pseudo nodes in free GIS software?

this might be the closest answer that i am looking for but this is related to lines and i need to delete them :(

examples of Pseudo Nodes in Polyline is illustrated below. (and the second picture, i need to apply this to polygons too)

an example of Pseudo Nodes in Polygon is illustrated below

so, the definition of Pseudo Nodes here is unnecessary nodes in a shape.
Also, with or without the pseudo nodes, the shape will still remain the same.

For instance, straight line only needs end nodes (2 nodes) but somehow there is a couple of extra nodes on the line (for some reason. maybe editting error), i would like to remove them.


texts: Using gdal_translate with a high number of gcp

texts:
I've a CSV file witch describes the coordinates of each pixels in a image, it is structured like this:

xcoord,ycoord,row,col
716657.6518375501,4784215.095395751,0,0
716657.6915582813,4784215.125764541,0,1
716657.7312790118,4784215.156133341,0,2
716657.7312790131,4784215.156133339,0,3
716657.810720475,4784215.216870922,0,4
....
....
716728.6399590928,4784196.171360771,1164,898
716728.6796798228,4784196.201729561,1164,899

Now I want to georeference the image using all of this points, I try to do the following script using gdal_translate:

def georeferencig (img,tab_coord):
    csvfile = open(tab_coord)
    tab = csv.DictReader (csvfile, delimiter=',')
    lista = []
    for row in tab:
        string_1 = '-gcp ' + repr(row['row']) + ' ' + repr(row['col']) + ' ' + repr(row['xcoord']) + ' ' + repr(row['ycoord']) + ' '
        lista.append (string_1)
    string = ''
    for part in lista:
        string = string + part
    output_name = img.split('.')[0] + '_refraction.tif'

translate = 'gdal_translate -of GTIFF -a_srs EPSG:32631 ' + string + img + ' ' + output_name
os.system(translate)

If I run this script the system gives me no results or errors, I try to use just a little part of the points writing for example:

def georeferencig (img,tab_coord):
    csvfile = open(tab_coord)
    tab = csv.DictReader (csvfile, delimiter=',')
    lista = []
    for row in tab:
        string_1 = '-gcp ' + repr(row['row']) + ' ' + repr(row['col']) + ' ' + repr(row['xcoord']) + ' ' + repr(row['ycoord']) + ' '
        lista.append (string_1)
    string = ''
    for part in lista:
        string = string + part
    output_name = img.split('.')[0] + '_refraction.tif'
    example = "-gcp '244' '495' '716683.4853030661' '4784220.244624769' -gcp '244' '496' '716683.5250237972' '4784220.274993559' -gcp '244' '497' '716683.5647445283' '4784220.305362349' -gcp '244' '498' '716683.6044652605' '4784220.335731139' -gcp '244' '499' '716683.6441859904' '4784220.36609993' -gcp '244' '500' '716683.6839067215' '4784220.3964687195' "

translate = 'gdal_translate -of GTIFF -a_srs EPSG:32631 ' + example + img + ' ' + output_name
os.system(translate)

And with only this points the script works fine.
Do you think that the problem is the high number of gcp written in 'string' or maybe some part of the script are not correct?
If I print 'string' I can see that is correct but has more or less 59.000.000 characters.


texts: Field Calculator case sensitive issue in QGIS graphical modeler

texts:
I use QGIS 3.10 with the graphical modeler and I'm facing an issue while creating new field in field calculator: it seems to be case insensitive thus it does not create the field I want.
Problem description follows:

Using Distance matrix algorithm, a layer is created with a field named Distance (and the name of the field can't be chosen)
I want to rename this field to distance to be compliant with further processing tools (REST API expecting a distance attribute)
I tried using Field Calculator to create a new field distance based on the value from the Distance field
When I run the model, the new field isn't created

If I choose a different name for the new field (like foo) everything runs fine. But how can I create 2 different fields with only case difference (then remove one of them).
I tried using Refactor fields as an alternative, but if I don't know the list of attributes in input layer, I can't properly use it.
Am I missing something in the configuration or is it an issue I should report on qgis github ?


texts: Prerequisites for running Visual Studio 2012 ArcGIS desktop application exe on any machine don't have ArcGIS installed

texts:
Scenario:

i created an exe application using Visual studio 2012 ArcGIS Map Control Application.
When running this application on any other machine which don't have ArcGIS installed on it. Its showing SystemIOFile exception and application don't open.
This application uses Esri.ArcGIS.Version and other dll's in it.

Any Solution?


texts: I have a point feature class (Address Points) and I am interested in creating a polygon feature class around groups of 500 address points. Is there a way I can do this programmatically, as opposed to manually creating grids?


texts:
I want to draw curves from a given point to another in Sample Flex Viewer like shown below.

Any help?


texts: Math formula for transforming from EPSG:4326 to EPSG:3857

texts:
Is there any math formula for transforming from EPSG:4326 to EPSG:3857 ?
i want to implement it in php code


texts: OL3 check constrain point drag to lineString?

texts:
I am wanting to allow the user to be able to drag a map marker, in OpenLayers 3, but constrain the dragging to the lineString that represents the route. Based on the code example on Stack Exchange I am able to move the marker and have an onChange handler, as follows:

geoMarker.on('change',function() {
    var point = scope.getGeometry().getCoordinates();
    if ( scope.map.getLayer('routeLayer').getSource().getFeaturesAtCoordinate(point) == 0 ) {
        // TODO if point is not in lineString, then prevent drag
    } else {
        // TODO get index of nearest point in routeLine and display
        //      location information
    }
});

I looked at the OpenLayers 3 examples, but nothing jumped out as how to go about doing this.


texts: Conceptually, what's the ideal way to develop detour maps for a large number of nodes

texts:
I have a road layer of roughly 30,000 links. On these links, I have a total of 15,000 nodes that are sorted in a point layer.

I'm trying to build a detour distance and routing around each asset in case that asset fails.

I have access to ArcGIS (not with network analyst) and QGIS latest build. I have access to PostGIS and pgRouting but I'd rather stay perhaps in network tools whether it's ArcPy, QGIS python or even straight python with shapely and netowkrX

Currently, I'm in the solution design stage. I'm curious as to how you'd recommend I build a "rerouting" engine that can efficiently output detour_distance and list_of_routes_of_detour

How would you recommend I build such a system?


texts: Connect a sub-model with multiple outputs to current model

texts:
I have a model in model builder that is clipping multiple feature classes to my study area. It is running with an iterator and creating 4 new clipped feature classes that are only represented by one output in model builder. Now I am trying to add this model as a sub model to run further tools on the individual outputs from my sub model. My problem is that the sub model only has one output and any tool that I attach to it will run for each feature class that was created from the sub model. Is there a way to have the sub model show each individual clipped feature class to connect them to separate tools?

For Clarity*

I'm currently running a bunch of clips to my study area at the start of my model that looks like this.

There area 4 feature classes being clipped in my data set. In an attempt to clean my model up, I created a sub model with an iterator to clip all those feature classes and it looks like this.

Now all my feature classes are run through the iterator, however, there is still only one output icon in the model. So I can not connect individual outputs from that iterator to the rest of the model. If I connect that final output, it will run the model for every layer that goes through the iterator.

The above model is showing how the sub model is represented in my current model, but only has one output icon and therefore I can't connect individual output features from the sub model to the current model.


texts: Automatically fill empty area between polygons if within certain distance using ArcMap?

texts: Suppose there are two polygons that are closer than 10m to each other. I would like to either combine those two polygons or create a new polygon that fills the empty area between them.


texts: Earth Engine Feature Collection to Image doesn't return values (can't work with Image)

texts:
I have a bunch of Polygons and build a Feature Collection out of them.

For testing the only property of each Feature is "Code" with the value 3:

import ee
import geopandas as gpd
ee.Initialize()
from geomet import wkt
from IPython.display import Image
import json

features = ee.List([])

x = 0
for i in range(0,len(gpd_frame)):
    g = wkt.loads(gpd_frame["geometry"][i].to_wkt())        #this part is for getting the geometry of the feature out of a dataframe
    fet = ee.Feature(g,{"Code":ee.Number(3)})
    features = features.add(fet)
    x = x+1
    if x > 500:                                                                #for testing
        break

col = ee.FeatureCollection(features)                              #FeatureCollection from list features

imgFC = col.reduceToImage(["Code"],"first")

imgFC.getInfo()

{'type': 'Image',
 'bands': [{'id': 'first',
   'data_type': {'type': 'PixelType',
    'precision': 'int',
    'min': -9223372036854775808,
    'max': 9223372036854775807},
   'crs': 'EPSG:4326',
   'crs_transform': [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]}]}

If I now want to look on the values with:

values = imgFC.select("first").reduceRegion(reducer=ee.Reducer.toList(),geometry=gebiet,maxPixels=1e8,scale=10).getInfo()

I get the error:

---------------------------------------------------------------------------
EEException                               Traceback (most recent call last)
<ipython-input-119-a69761c0332a> in <module>()
----> 1 values = imgFC.select("first").reduceRegion(reducer=ee.Reducer.toList(),geometry=gebiet,maxPixels=1e8,scale=10).getInfo()

~\Anaconda3\lib\site-packages\ee\computedobject.py in getInfo(self)
     93       The object can evaluate to anything.
     94     """
---> 95     return data.computeValue(self)
     96
     97   def encode(self, encoder):

~\Anaconda3\lib\site-packages\ee\data.py in computeValue(obj)
    488     The result of evaluating that object on the server.
    489   """
--> 490   return send_('/value', ({'json': obj.serialize(), 'json_format': 'v2'}))
    491
    492

~\Anaconda3\lib\site-packages\ee\data.py in send_(path, params, opt_method, opt_raw)
   1184       raise ee_exception.EEException('Invalid JSON: %s' % content)
   1185     if 'error' in json_content:
-> 1186       raise ee_exception.EEException(json_content['error']['message'])
   1187     if 'data' not in content:
   1188       raise ee_exception.EEException('Malformed response: ' + str(content))

EEException: An internal server error has occurred (ecbf2f1b7ee629dde4a2080f8856209a).

Has anyone an idea to solve this or another way to add features/polygons with values (later with different values for "Code") to another image?
As result the Image should have a band with the values of the features.


texts: DXF export from QGIS is not exporting all features

texts:
I'm trying to export multiple layers from QGIS 2.18 to DXF format for opening in CAD.

I am using the native DXF Exporter (Project -> DXF Export).

Here are the parameters I'm setting:

And the result in CAD does not contain all of this.

At a guess, it only seems to be exporting line layers and ignoring the polygons? Although this seems odd.

Can anyone suggest why I am losing features?


texts: I am trying to do a Euclidean distance analysis on my polyline shapefile. When I try to put the cell size as default when trying a Euclidean distance together with a given processing extent everything seems to work but the cells are too big to be process later with my other layers on weighted overlay. When I tried to put my cell size together with a given processing extent I get an error 99999. I already checked my projections and they are on the right WGS 1984 projection. I also tried giving the output path a different folder other than .gdb but it still doesn't work.


texts:
I would like to have a keymap in QGIS as in ArcGIS Desktop.
ArcGIS Desktop allows one to add two layers. You then load the main layer in the main map layout and another in the keymap and style them as preferred.

In QGIS, I load layers that I can group. Then on the layout, I create two maps (add two separate maps with different sizes). The main one and the key map. However, all layers loaded in the project shows up in the two maps added on the layout. How do I separate the two maps and style them accordingly. I tried to have two groups, style as prefered and then lock the key map, however after I reload the other layers, they still show on the locked key map.


texts:
I have a map with a few different layers. One of these layers is a 'cities' layer. I am trying to use selection->select by attribute-> Select: "Pop_CL"=8 or 9 or 10..... where Pop_CL is the population class (ten being cities with the largest population). I keep on getting error messages though!

Error 1: There was an error executing the query. An invalid SQL statement was used. I am selecting each # from "get unique values." No idea why I keep getting an error!


texts: Data Source of layer changing after copy/paste in Table of Contents of ArcMap?

texts:
I run into this situation all the time;

I often just copy layers and groups of layers in the TOC of one Map Document and paste then into the TOC of another Map Document.

For the most part this works fine but often enough the Data Source of layer changes when pasted. As a result the pasted layer has the red exclamation mark beside and no data is displayed until I manually set the Data Source.

I confirmed that the Data source of the layer I copy on one document changes when I paste the layer into another document.

This is something that has been happening since the early 8.x version and still an issue in version 10 SP3.  It's also not computer specific as i run into this issue on several computer.

Any ideas?

Layer after being pasted:

Source "layer being copied" Data Source:

Destination "after pasted in another document" Data Source:


texts:
I'm currently moving from using QGIS more as a visual tool to an analytical one. I am hoping to get better at optimal location distribution.
I have a dataset of census population data visualized as a shapefile. I'd like to optimally distribute buffer circles across the U.S. to maximize the capture of population within each buffer. So far, it looks like most plugins require a start point for a buffer and calculations are completed around that point. I'd like to flip that so that data determines which point to select.
How would I go about setting something like this up?


texts: Replace the existing data sources all MXD files with a new data source

texts:
I recently replaced my existing data server with a new server thus a new server name and path.
My existing data server is named skunk while my new server is tesclone.
I need to change all the datasources in my over 400 MXD files that currently point to skunk so that they point to tesclone.

The script has been modified using input and suggestions from members on this forum and now the script runs without error but unfortunately it still does not change the datasource.
I have tried using the print command to check if the scripts is able to access both the new datapath and the old datapath as well as the layer sources, all these were done successfully.

But still i do not understand why the line "mxd.findAndReplaceWorkspacePaths(oldPath,newDataPath,False)"  is not working

Modified script is shown below and the new generated output as well

#import the necessary python modules
import arcpy,os,csv

#setting folderpath
folderPath = r"./"

#setting environment variables
ws = arcpy.env.workspace = folderPath
arcpy.env.overwriteOutput = True

#Tell the script which data path needs to be updated
oldPath = r"\\skunk\GIS\Projects\Irene"

#Tell the script what the new data path is
newDataPath = r"\\isilonc1\dept\tnr\skunk\GIS\projects\Irene"

#### Stop editing the code here ####

#use arcpy to list all of the mxd files
def ListMXDSources(path,extension):
list_dir = []
CountList = []
MapList = []
list_dir = os.listdir(path)
count = 0
for paths, dirctory, files in os.walk(path):
for file in files:
    if file.endswith(extension): # eg: '.mxd'
      MapList.append(os.path.join(paths, file))
return(MapList)
print MapList
for m in MapList:
            count += 1
            mxd = arcpy.mapping.MapDocument(m)

#-------------Use below to call the function----------------#
mxdList = ListMXDSources(path=r"\\skunk\GIS\Projects\Irene",extension=".mxd")

#print the total number of mxd files
if len(mxdList) == 0:
 print "There are no mxd files found in this directory."
if len(mxdList) == 1:
 print "There is "+str(len(mxdList))+" mxd file found in this directory."
if len(mxdList) > 1:
 print "There are "+str(len(mxdList))+" mxd files found in this directory."

#set the count of the mxds to 0
count = 0
errorCount = 0

#starting the loop to check for each of the mxds
try:
for mxds in mxdList:
#set the mxd to the workspace with the mxd files
    mxd = arcpy.mapping.MapDocument(mxds)
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.supports("DATASOURCE"):
            theData = lyr.dataSource
            dataSourceCheck = lyr.dataSource[:-8]
        if dataSourceCheck == oldPath:
           mxd.findAndReplaceWorkspacePaths(oldPath,newDataPath,False)
           count += 1
           mxd.save()
        print mxds
    #del mxd
 except Exception as e:
print "The following error has occured:"
print e.message
print "Because there was an error, you should try re-running the script."
errorCount += 1

successCount = count-errorCount
print
print "The conversion process has finished and "+str(successCount)+" file(s) have been replaced."


texts: Aggregating or summing field by two or more criteria using QGIS?

texts:
I am using QGIS 2.18.

I'm trying to sum fields based on two criteria with filed calculator.

I have the "ID" field and the "ZONE" field (ID is secondary key, repited for other table). ZONE field has a value 01;02;03.I want to summarize the field "VALUE"  based on the ID field, and ZONE field.

aggregate( 'ID', 'sum' , "VALUE" , "ZONE" = 01)

aggregate( 'ID', 'sum' , "VALUE" , "ZONE" = 02) ect.


texts: Finding the best image per year of Landsat image collection in Google Earth Engine

texts:
I am trying to create a series of images with the best image per year using Landsat 7 from 2000-2017 around the legal Amazon (using RBG bands 1-3). The code I have so far pulls up over 70,000 images over the period. It seems like I should be looping through images in each year and finding the least cloudy image, but I know you're not supposed to write loops in GGE.

I'm really new to this platform so I may be missing something obvious!

The end goal is to create a timelapse that I will overlay with my own dataset. The code below is what I have so far, and I have no idea how to proceed from here:

// create polygon around legal amazon
var Amaz = ee.Geometry.Polygon([[-43.0, -18.0], [-69.0, -18.0], [-69.0,5.0], [-43.0, 5.0]]);
Map.addLayer(Amaz)

//dates of interest
var start = ee.Date('2000-01-01');
var end = ee.Date('2018-01-01');

// create image collection

var Amazon_images = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
.filterBounds(Amaz)
.filterDate(start, end)
.sort('CLOUD_COVER', false);

// Get number of images
var count = Amazon_images.size();
print('size of collection Amazon_images', count);


texts: Merging two atributes columns in ArcGIS for Desktop?

texts:
I have a shapefile that is a result of a merge between a land use chart (with two gaps) and two files to fill the gaps, but the atributes in the Land Use Chart are designated by "LEG" and in the Gap Charts are designated by "COS".
In the merged chart, I got two columns that I would like to merge.
They are both listed as "Text".
I already tried the Field Calculator but the SQL code was invalid.
Could anyone help me?

I am using ArcGIS 10.1 for Desktop.


texts:
When using the term "feature" in GIS, what is often (if not always) meant is a geographical feature which "is a representation of real world phenomenon associated with a location relative to the Earth" (ISO 19125-2:2004). Are there non-geographical features and if so, are there any examples of such?

For example, could a compass or a scale bar be considered non-geographical features?


texts: Closest feature at angle using ArcGIS for Desktop?

texts:
I have points representing where street light poles meet the ground. I also have the angle which the mast arm is facing. We use this to symbolize the features by showing the mast arm over the appropriate road- sometimes it is not perpendicular to the closest road. I need to associate attributes in the street centerline with the point.

The near tool doesn't do it because you can't specify a specific angle.


texts: How can I add my map to my personal blog?

texts: I have been able to deploy the map on the web, so I can access the html file and all the javascript together, but can I put it on my personal website? My blog has been deployed using Ruby on Rails


texts: How to randomly select point within grid (QGiS)

texts:
I have a series of occurrence points that I am trying to downsample to 1/20 km to avoid sampling bias. So far, I have made a grid of ~20km squares and have filtered them down to squares with more than one occurrence point in them.
I want to figure out how to randomly select a single point within each grid square and delete the rest in the square. I tried joining the point vector and grid vector layers as was suggested in this post, but whenever I try to filter by longitude or ID, it selects every point in every square except one.
Does anyone know either what I'm doing wrong that's preventing me from selecting points randomly or a simpler method to doing this?


texts: Geotools : A layer to draw the traffic volume of a road

texts:
I have developed a standalone swing application loaded with a shapefile of a road network using geotools (v.14.3). Now I would like to implement a feature where the user may zoom the map, and at a pre-set zoom level, the traffic volumes of the currently visible road would be drawn close to the road (Here's a sample, instead of letters, there would be the traffic volume in a number). Traffic volume is a number which the user would need to take decisions.

I've seen an example here where they discuss about drawing circles and other polygons on a separate layer on a JmapPane.

How can I draw the traffic volume of the currently zoomed road?


texts:
I have a 3D multipatch feature class in a file geodatabase using ArcScene 10.4.1. I would like to export this 3D data to CAD using "Export to CAD (Conversion)".

The tools runs, but the resulting DWG-file is only 12kb and is empty.

Is it even possible to export 3D multipatch-files to CAD using ArcGIS?

Just exporting the 2D footprints to CAD is working.

The geometry is fine.


texts: Adding new field values to layer from another layer

texts:
I have two layers that share a Unique "ID" field. Can I easily update the field values in one layer with the field values in the other layer using the Unique "ID" as the link field? As the layers also sit on top of each other I suppose I could alternatively use the 'equals' operator in 'Join Attributes by location' but I thought using a fixed field value would be more reliable. One way seems to be to join the two layers and then to perform a field update but I thought it might be possible in one step with the right function.
I am aware that I can join the layers, but I need to get the data into the original layer itself so they can be handled separately.


texts: How to use QGIS to create GeoTIFF contour map to import into Mapbox?

texts:
The final product I'm trying to achieve is a Mapbox interactive map with a colored contour layer.

I have a CSV file that contains rows with "LATITUDE,LONGITUDE,VALUE" in a 0.01-degree grid. I want to create a contour map by interpolating the points and creating a raster map. From my research, I have determined that QGIS should be able to do this, but I'm not quite sure how.

I believe I need to output a GeoTIFF to then import into Mapbox. I found this potential solution: https://stackoverflow.com/questions/26715804/how-to-create-interactive-online-map-with-contour-filled-plot-layer, but I'm hoping to get more information on Steps 1 and 2 of that solution.


texts: Creating a new field having the same type as selected field?

texts:
I set up a QgsFieldComboBox for a Plugin. I know I can read the selected fieldname the following:
self.dlg.MyQgsFieldComboBox.currentField()

Studying the docs, I cannot find an option to read the currentFieldType() directly. So I know I can get field types of my selected layer the following:
self.mylayerfieldtypes = [field.typeName() for field in self.myselectedlayer.fields()]

I could stitch this together the following, to get the field type only of my selected field:
myfieldname = self.dlg.MyQgsFieldComboBox.currentField()

for field in self.myselectedlayer.fields():
   if field.name() == myfieldname:
       myfieldtype = field.typeName()

But my issue now is, when creating a new layer having a field with the same type as my selected field, I cannot just use myfieldtype for it like
QgsField("MyNewFieldHavingTheSameTypeAsMyOldField", myfieldtype)

How can I create a new field, having the same type as my selected field? Copying the field to my new layer is not an option.


texts: Create extension PostGIS throws unknown error 193

texts:
I have execute all actions like this
How to install PostGIS on Windows?.
But when I create extension postgis;

I have got an error message :

  "PostgreSQL/9.3/lib/postgis-2.1.dll » : unknown error 193".

How can I resolve this?


texts:
I'm running GeoServer 2.15.2 on Ubuntu 18.04.3 with the latest ubuntugis-unstable release packages (GDAL 3.0.2, PostGIS 3.0.0) and cannot get the Image Mosaic *.jar to load. I followed the instructions of https://docs.geoserver.org/stable/en/user/data/raster/imagemosaicjdbc.html

The pertinent files are loaded into /usr/local/geonode/geoserver/geoserver/WEB-INF/lib

gt-imagemosaic-jdbc-21.2.jar
gt-imagepyramid-21.2.jar

The pertinent xml configurations are in /var/lib/geoserver_data/coverages and are

connect.pgraster.xml.inc  ned_13.pgraster.xml  ned_1m.pgraster.xml
mapping.pgraster.xml.inc  ned_19.pgraster.xml  ned_1.pgraster.xml

The contents of connect.pgraster.xml.inc are

<connect>
   <dstype value="DBCP"/>
   <username value="geonode"/>
   <password value="password"/>
   <jdbcUrl value="jdbc:postgresql://localhost:5432/spl_gis"/>
   <driverClassName value="org.postgresql.Driver"/>
   <maxActive value="10"/>
   <maxIdle value="0"/>
</connect>

The contents of mapping.pgraster.xml.inc are

<spatialExtension name="pgraster"/>
<mapping>
    <masterTable name="MOSAIC" >
        <coverageNameAttribute name="NAME"/>
        <maxXAttribute name="maxX"/>
        <maxYAttribute name="maxY"/>
        <minXAttribute name="minX"/>
        <minYAttribute name="minY"/>
        <resXAttribute name="resX"/>
        <resYAttribute name="resY"/>
        <tileTableNameAtribute  name="TileTable" />
    </masterTable>
    <tileTable>
        <blobAttributeName name="rast" />
    </tileTable>
</mapping>

The contents of ned_1.pgraster.xml are

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ImageMosaicJDBCConfig [
    <!ENTITY mapping PUBLIC "mapping"  "mapping.pgraster.xml.inc">
    <!ENTITY connect PUBLIC "connect"  "connect.pgraster.xml.inc">
]>

<config version="1.0">
    <coverageName name="ned_1"/>
    <coordsys name="EPSG:4269"/>
    <!-- interpolation 1 = nearest neighbor, 2 = bipolar, 3 = bicubic -->
    <scaleop  interpolation="1"/>
    <axisOrder ignore="false"/>
          &mapping;
          &connect;
</config>

The above files should only be needed once the ImageMosaicJDBC is selected. As such, these are untested...

Once tomcat has been restarted the two *.jar modules do not appear in the Server Status - Modules list neither does the ImageMosaicJDBC option appear when adding a new data source.

GeoServer accesses the distro PostgreSQL 12 just fine for all of the vector data. So, I'm not sure what/which config I'm missing. Any pointers?

EDIT
There are a number of warnings during startup that aren't causing any known issues. There was a single warning returned for image mosaic:

WARN [gce.imagemosaic] - Unable to set ordering between tiff readers spi


texts:
>>> import arcpy
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Python27\lib\arcpy\__init__.py", line 21, in <module>
    from arcpy.geoprocessing import gp
  File "C:\Python27\lib\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from _base import *
  File "C:\Python27\lib\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting
ImportError: DLL load failed: The specified module could not be found.

This is the error I get when i try to import arcpy. Any suggestions?


texts: How to enable Arabic/UTF8 labelling in OpenLayers?

texts:
I have data in GeoJSON that contains Arabic in some attribute when reading the attribute

info.push(features[i].get('Name'));

I get this:

Unreadable how to enable UTF8 or I don't know...


texts: Raster processing with Python or Arc raster calculator

texts:
I have 4-band NAIP imagery that require the following processing steps:

Convert to 8 bit unsigned B/W
Reverse (note not inverse) the image, so that dark values appear
light and light values appear dark.  Moderate values should remain
relatively unchanged.

Could someone please point me in the right direction for a Python or Arc raster calculator approach?


texts: Project from GCS_WGS_1984 to WGS_1984_Web_Mercator_Auxiliary_Sphere

texts:
What is the default transformation method to project data from from GCS_WGS_1984(WIKI 4326) to WGS_1984_Web_Mercator_Auxiliary_Sphere(WIKI 3857)?

ArcMap 10.3.1


texts: Rules for merging map markers at scale in QGIS

texts:
I have a map of Australia, plotting cafes.  As I change the scale, these markers cluster too tightly to be useful.  Point displacement is also not useful, and a heat map is too indistinct.  What I'd like to do have those markers collapse into one marker, which displays a count of nearby markers.  For example, I'd like, when the zoom reaches 5000:1, for a single map marker to display "78" - the number of individual markers that are close to each other at that scale at Sydney's CBD.

I expect this is going to involve a complex render rule.  Is this the best way to do this?  If so can you give me a pointer as to how to call the scale factor in the rule, and how to set the point to which the markers collapse?


texts: Getting bounding box of a selected region through PyQGIS?

texts:
Assume I have selected several features of a layer in QGIS and I would like to get the bounding information of the whole selected area through python console like the bounding of yellow regions the figure shows.

What is the easiest method to do it?

p.s. I have tried to use layer.selectedFeatures() but it just extract the list of selected features but cannot directly get the bounding info.