texts: Network Analysis Shortest path with original line segments as output

texts:
I've managed to batch run the shortest path algorithm and the shortest route (point to layer) in the output is perfect, except that I really need the individual attribute tables of the line segments in the output. Any idea if there is a way to produce something like this?

I'm basically looking to create a new version of my line input layer with the all original segments/attributes intact, but with the ID of whichever start point (x however many, if the same line feature is used by more than one start/end point combination) was used for that output line feature.

Right now, I’m getting complete line from a location to z location, but I need to know which ‘roads’ make up that full path.

Edit: added screenshots of data.


texts: MODIS calculations with non cloud contaminated pixels

texts:
I have 11 files showing 8day daytime LST temps in .tif format from MODIS.

Is there any way to calculate the average 8day daytime temperature for march 2013 for example but using in the calculations only pixels which are not cloud-contaminated in any of the 11 files that participate in this process?


texts: Copying data from one attribute field to another using ArcGIS Field Calculator?

texts:
I am standardizing several shapefiles.

How do I go about copying data from one attribute field to a new attribute field i.e. would I use the calculate field tool?


texts: How to change the Fieldname(s) in the attribute table

texts:
I joined some data from excel to arcmap, but the added fieldnames from the excel table (land, Germany, field3, field4, ...) aren't the fieldnames which I want to use. For example land=district and Germany has to be 2006. When I only change the alias name, the name removes when I turn my PC off... Is there a way to change them permanently.

I found out that the data was visable at the attribute table, but it was not at the fields in the shapefile properties. By exporting the joined shapefile to a new shapefile this problem is fixed, but the fieldname(s) are still not editable.

I saw you can edit the names at the shapefile properties, but when I change them it gives a "Alter field name Failed" error.

Does somebody know a solution?


texts: Georeferencing an image in R. When I change the extent, x and y resolutions become different

texts:
I am georeferencing a number of images using R. I can read the image in fine, and set an extent:

library(raster)
library(rgdal)
library(maptools)

## Read in the image as a raster
my.raster <- raster("MyImage.tif")

## Read in previously gecoded maps
my.polygons<-readOGR("Poly.shp")

## Fix the extent
extent(my.raster)<-extent(my.polygons)
projection(my.raster) <- "+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0"

But when I look at the file I see that the x and y resolutions are off:

class       : RasterLayer
band        : 1  (of  3  bands)
dimensions  : 1935, 2511, 4858785  (nrow, ncol, ncell)
resolution  : 0.0006940597, 0.0003798972  (x, y)
extent      : -91.81082, -90.06803, 13.84329, 14.57839  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0

and then when I start manually changing the extent to match the shapefile I have, it only gets worse. If I resample(from the raster package), it gets worse again.

Can I georeference an image and keep the resolutions consistent so that when I export it to ArcGIS or whatever, it is properly placed?

I am not providing files, because the project has confidentiality issues.


texts: Polygon splitting and merging to another one without overlapping using ArcGIS Desktop?

texts:
I have a layer of different types of polygons. Some are fields, other are roads. Sometime a road polygon is entering a field polygon. As presented on the capture, the yellow road (polygon 1) is getting inside the blue field (polygon 2).

I am trying to split the road in 2 pieces 1A and 1B, at their edges intersection (following the red line), so I could then merge and dissolve the part 1A of the road to the field (polygon 2) it is contained within.

The goal is to have a layer with only fields polygons. Roads should be erased, and all gaps filled and merged to the largest adjacent polygon.

As there are no overlapping, I was looking of a way to cut/split the road polygon, before merging it to the fields. But as there is a very important amount of data, this process should be done automatically, and not manually by identifying visually edges and splitting using polylines.

Is there a way you can think of, using ArcGIS tools (perhaps production mapping?) or any other tool or method that I could use?


texts:
I have a geodatabase that I am adding to a QGIS project.  I read the documentation (below) but I don't understand what parameters I use for "a.the_geom" and "b.geom".   Is there a way to clip with a bounding box?  Or, what is the best way to reduce my huge data to just the small portion that I need?

  16.5.9. Clipping
  To clip a subpart of your data you can use the ST_INTERSECT() function. To avoid empty geometries, use:

  where not st_isempty(st_intersection(a.the_geom, b.the_geom))

I am using QGIS 2.8.2-Wien, a portable version. I'm running QGIS from an external hard drive connected via USB.I've been opening a new project and then using Layer > Add Layer > Add Vector Layer > Source Type File -- and adding some of the files files that are in a .gdb geodatabase file.

It sounds like I can use the Select Features by Area or Single Click tool to select a rectangle and then I should be able to right-click or somehow save the selected area as a shapefile. However, when I select an area, as soon as I let go of the mouse button, I get the crash dump shown in the image below and QGIS shuts down. I will try using an older more stable version of QGIS.

The method I ended up using is probably only one of many ways to cut/crop/clip a small portion of data from a large data set.
The first thing I did was switch to an older, stable, version of portable QGIS - 2.4.0-Chugiak.
In the Layers Panel, I right-clicked on the layer, clicked Save-As...   a "Save vector layer as..." window popped up which let me choose which format to save as from a drop-down menu at the top and there is a checkbox for Extent which allowed me to put in bounding box coordinates.  I ended up with a shapefile with only the area that I wanted.


texts: Geometry has too many edges error on GEE

texts:
I was trying to perform an intersection geometry operation on an asset I uploaded which I converted to a feature collection and I keep getting the "intersection: Layer error: Collection.geometry: Geometry has too many edges (18183861 > 2000000)." error regardless of size of polygon.
Link to my work
https://code.earthengine.google.com/73f647cc64d430d60b69a021127c148f
// converting the asset to feature collection

var irrigation = ee.FeatureCollection(Irrigation);

// Simplify each feature in the collection, by mapping the
// .simplify() function over it.
//var simplifiedirrigation = irrigation.map(function(irrigation) {
//  return irrigation.simplify({maxError: 100});
//});

// create a polygon

var polygon = ee.Geometry.Polygon([
  [[5.364373355349645,7.161499312988462], [5.364373355349645,7.161499312988462], [5.891717105349645,7.597311534607618], [5.891717105349645,7.597311534607618], [5.364373355349645,7.161499312988462]]
]);

// Compute the intersection.

var intersection = polygon.intersection(Irrigation);

// displaying data

Map.addLayer(irrigation, {}, 'Irrigation Potential');
Map.addLayer(intersection, {}, 'intersection');
Map.setCenter(5.7621,7.1942, 8);


texts: Cutting polygons when it becomes too narrow

texts:
I'm using QGIS 3.4.5.
I'm working with a polygon layer that has a large set of polygons. Sometimes the polygons have a, for their purpose, irrational segmentation, with long and narrow parts. I would like to find a way to cut the polygon when it's becoming "too narrow" (green lines in the image). This would then create two polygons, of which the smaller polygon can be filtered out by area.

I tried to convert it to lines and then run the plugin "Check geometry" to search for a maximum angle degree, but that did only remove the very tip of the narrow part. I also did try this with overlapping areas but couldn't get a good result.
As an alternative to the more ideal approach above, I've also tried to create a negative buffer to eliminate the narrow parts, and then do a regular +buffer to achieve a polygon that can be cut with the original layer. This procedure is OK, but will result on overlapping polygons, which I'm not sure how to handle correctly. When I run v.clean to create individual polygons of the overlapping areas they sometimes have similar or higher area as the correct polygons, why they can't be eliminated by filtering for area.
Is there any way/method to either:

Cut narrow parts of a polygon at a certain narrowness/distance, or
Create a mean common polygon border between two overlapping polygons such as this green dotted line:

Any ideas?


texts:
I have thousands of vectors (out of several million) which I need to flip, for directionally purposes (Network Analyst).

The vectors in question are annotated that their directionality has been reversed.

Anyone have or know of a way to flip multiple vectors at once or have any code which would do the same?


texts:
How do I filter data by date. If I want to filter, I can see the old filtering too. The layers are stacked on top of each other. How do I get the current date range.

<div id="time-ranges" class="leaflet-bar">
    <form action="#" id="form">
      <label for="from">
        From
        <input id="from" type="date" value="1998-07-01" name="from">
      </label>
      <label for="to">
        To
        <input  id="to" type="date" value= "2019-12-19" name="to">
      </label>
      <input id = "sub" type="submit" value="Update">
    </form>
  </div>
$("#form").submit(function(e){
        e.preventDefault();
        var inf = L.geoJSON(inf_st,{
            filter: function (feature, layer) {

                return new Date(feature.properties.fest_datum) >= new Date(document.getElementById("from").value) && new Date(feature.properties.fest_datum) <= new Date(document.getElementById("to").value);

            }

        }).addTo(map);

    });


texts: Symbolizing using RGB values from Excel in ArcGIS Desktop?

texts:
I have an excel file with 2 attribute columns, 1-feature info 2-Number of color (RGB) that the current feature should be symbolize. This file file can be joined with the layer using the "feature info" field which is unique (join by attribute table). So i have as attribute the color number attribute for every polygon feature.

How can  i use this this information to symbolize every feature with their attributed color already imported from the excel sheet?


texts: Create Geometric Network Geodatabase Failed

texts:
I recently attempted to run the Create Geometric Network Geodatabase tool which is part of a water outage template download and process offered via the ArcGIS Solutions Deployment Tool within ArcGIS Pro. The tool will run almost to completion but then fail immediately at the end.

Clearly this is missing from the AcrPy module but I'm confused as to why?

The following is my error message:

Traceback (most recent call last):
  File "<string>", line 235, in execute
  File "<string>", line 293, in run
  File "Y:\SYSTEM_SUPPORT\PILOT\Water_Outage_Template\Water_Outage_Template\WaterOutageComponents\Tools\ToolData\gnsolutions\publish.py", line 237, in main
    gn = self.create_network()
  File "Y:\SYSTEM_SUPPORT\PILOT\Water_Outage_Template\Water_Outage_Template\WaterOutageComponents\Tools\ToolData\gnsolutions\publish.py", line 99, in create_network
    result = arcpy.CreateGeometricNetwork_management(in_feature_dataset=self.dataset,
AttributeError: module 'arcpy' has no attribute 'CreateGeometricNetwork_management'
 Failed to execute (CreateGeometricNetworkGDB).


texts:
I have a query running as below:

UPDATE datasets.riverside_houses
SET easting = addressbase.x_coordinate
FROM datasets.addressbase
WHERE addressbase.building_number = riverside_houses.house_number OR addressbase.building_name = riverside_houses.house_number AND addressbase.thoroughfare = riverside_houses.address_line_1 AND addressbase.postcode = riverside_houses.postcode;

UPDATE datasets.riverside_houses
SET northing = addressbase.y_coordinate
FROM datasets.addressbase
WHERE addressbase.building_number = riverside_houses.house_number OR addressbase.building_name = riverside_houses.house_number AND addressbase.thoroughfare = riverside_houses.address_line_1 AND addressbase.postcode = riverside_houses.postcode

The table 'addressbase' is huge (it estimates over 30 million rows).

Essentially I have a table of addresses (riverside_houses - split into separate fields, same as addressbase) that I need to match to addressbase and return the easting/northing. I require it to check multiple fields as the condition, otherwise it could pull through all kinds of irrelevant results.

The code above comes from some research on stack, adapted to meet the above requirement. My small scale testing on sample data worked perfectly. This is my first big test and it seems to be struggling.

Is there a more efficient way to perform this action or do I just need to accept it's a big table and will take a long time?


texts:
I am somewhat confused on the meaning of the Overlap flag which is present in las spec 1.4.
The ASPRS LAS Version 1.4 Spec states:

New to LAS 1.4 is the fourth flag, OVERLAP. It is frequently desired to thin data in the overlap region of lidar flight lines (or, in the case of mobile, "takes"). This operation is used to either improve the statistical accuracy of the lidar data or to simply thin the multiple density data in these regions. Using a Class to indicate overlap (the technique prior to LAS 1.4) caused problems in that one loses the actual classification of the point (e.g.,Ground). This leads to problems in analyzing accuracy in the overlap region. What exactly comprises points that should be designated as Overlap will be determined (as it is now) by the procurer of the lidar data. For example, USGS has an extensive description of how points should be designated as Overlap.

The Lastools lasoverage readme states:

Below an explaination based on what Karl Heidemann of the USGS once told me regarding "overlap" versus "overage" points:
"Overlap" applies to any point in one flightline or swath that is anywhere within the boundary of another flightline. Whether or not a point is overlap is determined, basically, at the time of flight. It is a fundamentally immutable characteristic.
"Overage" applies to any point not in the "tenderloin" of the flightline or swath, that is, the central core of the flightline that - when combined with the others - creates a non-overlapped and gapless coverage of the surface.  Whether or not a point is overage is determined by the operator (or software) that defines the single-coverage tenderloins of the flightlines. Hence, it is a subjective characteristic -- it depends solely on where somebody (here: lasoverage.exe) chooses to draw the boundary.
Unfortunately the LAS specification uses the word "overlap" where it should be using the word "overage".

It seems the LAS Spec is using the overlap flag to refer to "overage" as defined in the lastools readme since the "overlap" points are subjectively defined by the "procurer" of the lidar data.  So, how should overlap points be handled? Should they be removed or used?


texts: To compare and merge two point layers based on lat long

texts:
I have 2 shapefiles of point type.  Both layers have lat/long associated with each point. In one shapefile, I have about 23 hundred thousand points and in another about 20k points.
Because the first one has 23 hundred thousand points or features, it is very difficult to process the data or compute any algorithms on it.

Hence I want to somehow compare both shapefiles such that i get a layer with only those points which are in common in both the shapefiles.
The shapefiles just have one common attribute that is geom (lat long) between them.

I am using WGS 4326 CRS .


texts: How to use Google Maps for more than just an overlay

texts:
I'm using ArcGIS Desktop 10.1 SP1 and I have defined the annual risk of accidents for each segment of a city road network.

I would like to know if anyone tried to do more than overlay on Google Maps or similar...

For instance, is there a way of using their services to return a path from A to B, based on my data? We could call it, "return the safest path from A to B", instead of the traditional shortest/fastest one?

I don't know if it's relevant or not, but my network doesn't completely match OSM data in ArcGIS (it was created that way and with a different datum ... so, it may be also due to the geographic transformation applied).

If there is a way, can it be achieved with my current license?


texts: Modifying CSV file in QGIS before joining it to shapefile?

texts:
I've imported a Delimited Text File that I want to join to a shapefile.

Unfortunately, the join field has a ' character before the actual string.

I wanted to delete that character, but I can't edit the file because I get:

  This layer does not support adding new provider fields

Consequently, I cannot edit the file (and yes, I toggled the layer editing icon on and tried editing it after saving the imported CSV with a new name).

Is there a way to get around this?


texts: How to programmatically distinguish a Sentinel2 partial and full tile?

texts:
I am downloading Sentinel-2 (and Landsat) images from the google public cloud using this script: https://github.com/vascobnunes/fetchLandsatSentinelFromGoogleCloud

This works fine for me, but in some dates I'm getting images like this:

which is not good for what I want.

I tried to find some info in the metadata where I could distinguish between a partial and a full tile but I am not finding any.


texts: Mapserver / QGis and DEM GeoTIFF different nodata values behaviour

texts:
I Have a problem with novalue datas in my DEM… This is the workflow I did:

my DEM datas are in multiple Arcgis ASCII grids… Each .asc files have a header containing "NODATA_value  -99999.00". I opened the in QGIS 3.0.3 and merged them in QGIS
I save my new layer in GeoTiff format, close QGIS and reopen my new GeoTIFF. All the "nodata" values have a value of -99999.
I process my nodata values with gdal_translate -co TILED=YES -of GTiff -a_nodata "-99999" input.tif outputNaN.tif
On QGIS the outputNaN.tif has a correct behaviour, nodata are "no data"
I created a Mapserver 7.0.7 map file:

MAP
  OUTPUTFORMAT
      NAME "GTiff"
      DRIVER GDAL/GTiff
      MIMETYPE "image/tiff"
      IMAGEMODE INT16
      EXTENSION "tif"
  END

  NAME "Modèle numérique de terrain"
    EXTENT  -6.0001389 39.9998611 10.0001389 52.0001389
    MAXSIZE 4096
    PROJECTION
        "init=epsg:2154"
  END

  WEB
    METADATA
        "wcs_label"           "WCS Server"
        "wcs_onlineresource"  "http://localhost/cgi-bin/mapserv?map=/tmp/testNaN.map&"
        "wcs_enable_request" "*"
            "wcs_fees"            "none"
            "wcs_accessconstraints"    "none"
            "wcs_keywordlist"          "wcs,dem"
            "wcs_metadatalink_type"    "TC211"
            "wcs_metadatalink_format"  "text/plain"
            "wcs_metadatalink_href"    "http://localhost"
            "wcs_address"              ""
            "wcs_city"                 "Paris"
            "wcs_stateorprovince"      "IDF"
            "wcs_postcode"             "75000"
            "wcs_country"              "France"
            "wcs_contactelectronicmailaddress" "admin@gnone.world"
            "wcs_contactperson"            "me"
            "wcs_contactorganization"      "unemployed"
            "wcs_contactposition"          "manager"
            "wcs_contactvoicetelephone"    ""
            "wcs_contactfacimiletelephone" ""
            "wms_title"           "WMS Server"
            "wms_onlineresource"  "http://localhost/cgi-bin/mapserv?map=/tmp/testNaN.map&"
            "wms_srs"             "EPSG:2154"
            "wms_enable_request"  "*"
            "wfs_title"          "WFS Server"
            "wfs_onlineresource" "http://localhost/cgi-bin/mapserv?map=/tmp/testNaN.map&" ## Recommended
            "wfs_srs"            "EPSG:2154 EPSG:4326 EPSG:5698" ## Recommended
            "wfs_abstract"       "WFS Server"
            "wfs_enable_request" "*"
      END
    END

  LAYER
        NAME TESTNAN
          DATA "/tmp/outputNaN.tif"
        STATUS ON
        TYPE RASTER
        PROCESSING "SCALE=AUTO"
        PROJECTION
           "init=epsg:2154"
        END
    METADATA
      "wms_title"         "DEM 75m"
      "wcs_description" "Modèle numérique de terrain"
      "wcs_name" "DEM 75m"
      "wcs_label" "DEM (rés: 75m)"
      "wcs_srs" "EPSG:2154"
      "wcs_enable_request" "*"
      "wcs_bandcount" "1"
      "wcs_band_names" "Altitude"
      "Altitude_band_description" "Altitude en mètres"
    END
  END
END

But this layer has now -32768 value where a "nodata" must be found. I tried to play with PROCESSING "NODATA=OFF" or PROCESSING "NODATA=-32768"but without success.

I know how to set in the layer properties in QGIS a transparency of 100% where pixel=-32768 . But I'd prefer to find a better solution. How can achieve this?


texts: Converting a GeoTIFF image from RGB to gray scale using Rasterio or GDAL in Python

texts:
I was using PIL to convert some images from RGB to gray scale in Python, but this code does not work well for GeoTIFF files. It seems that PIL does not seem to work with with GeoTIFF files in particular.
I was using the code:
from PIL import Image

t = Image.open('testfile.tif')
gray_image = image.convert('L')

So I was wondering if there is a way to load an image in say rasterio and then convert it to a  gray scale array. Now I could manually do this by writing a simple function to matrix multiply the image array, BUT the files are pretty large and so I run into memory issues. So I would likely need to do a windowed read on the file.
If there is a convenience function in rasterio or GDAL for this, could someone point it out. I searched but could not find anything. Otherwise I can try to work on creating a windowed read.


texts: ArcGIS Field Calculator nested elif is None not working

texts:
I am calculating a field based on two other fields in ArcMap 10.2 using python.

My Pre-Logic Script Code is

  def myfunc(legend,Btype):

if(legend == 'Apple'):
    return 'Banana'
elif(legend == 'Cranberry'):
    return 'Date'
elif(legend == 'Elderberry'):
    if(Btype.startswith('Fig')):
        return 'Grape'
    elif(Btype.startswith('Lychee')):
        return 'Mango'
    elif Btype is None:
        return 'Olive'
    else:
        return 'Pear'
else:
    return 'Raspberry'

Mostly this function works, however in the case that legend is "Elderberry" and Btype is "Null" the function returns "Null" rather than "Olive"

The problem only occurs when I nest the "is None" in the if statement.


texts: Replacing values in multiple columns using QGIS Field Calculator?

texts:
In QGIS I want to search for a specific case and if this is true then replace the text in multiple columns. Based on Use updatcursor to replace values in multiple fields, Update multiple fields/columns in QGIS and Elseif Conditional Statement in QGIS Field Calculator

I have

# definition of the function which calculates new values
def upd_value(sfield, svalue, outvalue):
    # criteria
    if search<>svalue:
     outvalue=svalue
     return outvalue

# get the active layer
layer = iface.activeLayer()
prov = layer.dataProvider()

# get the fields
fnm = prov.fieldNameMap()

# get an iterator for the features
feats = prov.getFeatures()

# calculate new values and update fields with fieldnames

I am not sure how to end this. Basically I need the following

When
 lower("Given_name")<>lower( "z_Given Names") and 'Surname' is not Null
Then
 "Given_Name" field should be filled with "z_Given_Names"
 "Surname" field should be filled with "z_Surname"
(lots of other fields which I can add in once I know what to do)
 "Status" with 'Updated fields from z table as names not same'

When
 "Status" = 'V' and "Surname" is Null
Then
 "Given_Name" field should be filled with "z_Given_Names"
 "Surname" field should be filled with  "z_Surname"
 "Status" with 'New'


texts: Using ArcPy with Python 3.7 instead of default Python 2.7 that comes with ArcMap 10.6.1?

texts:
I want to use ArcPy with Python 3.7. I tried setting the paths to arcpy but Python 3.7 doesn't import arcpy properly.
Is it possible to use a Python version outside of default version that ships with ArcMap 10.6?


texts: Avoid CORS problems using javascript code loading layers in OpenLayers

texts:
I am putting together a basic proof of concept for a client using OpenLayers to compare against Esri Javascript.

I understand that the medium to a long-term solution is to enable access on the server I am trying to gain data from. At this stage I am just prototyping and so wanted to load in a layer using the code, I probably will be using a proxy longer term anyway.

My starter for ten code is as follows:

<!DOCTYPE html>

    <title>Tiled ArcGIS MapServer</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v5.1.3/css/ol.css" type="text/css">

Presumably, I can add something from within the Vector source class to either use a proxy, or ignore any CORS issues?

The site is hosted in IIS BTW.

The error is :

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource
  here.
  (Reason: CORS header ‘Access-Control-Allow-Origin’ missing)

I have to admit I'm finding it difficult to find what I need in the Openlayers docs, couldn't find any demos using a proxy.

I have also found lots of posts around enabling CORS on the target server, this is not what I want to do right now at a sandbox prototyping stage. I am yet to persuade the client to use openlayers over Esri javascript (We will be wanting to manly use KML files from a server which we will NOT want to make public facing which is an issue with the ESRI control as it bounces via ArcGIS.com for some completely hard to grasp reason), so just want a quick and dirty prototype.


texts:
Edit:

Solved it. When I was creating my new zip/kmz my code was creating various folders within the zip/kmz.

Changed:
NewZip.write(addedfile)

to:
NewZip.write(addedfile, newfl)

/Edit

I have KMZs that need changes to their points. The changes are simply updating their attributes - changing the values in the tables that appear when a point is clicked.

My method involves extracting the files from the KMZ to a temp location, making a copy of the "doc.kml", digging through the copied "doc.kml" for instances of text that needs to be changed, and then renaming the old KMZ, creating a new KMZ with the same name as the original, and adding the all the files (including my edited "doc.kml") to the new KMZ.

I have automated this process with a pthython script, but when I run it, any KMZs that gets updated have broken point icons. i can't find anything that would cause this - my changes are quite minimal and do not replace and obvious KML. Changes are formatted: "PR_A_1567_20" to "PR_AB_1567_20", so essentially I am adding an extra letter into a few lines of the KML.

Anyone know why this might break icons? Everything else still functions perfectly.

Thanks.

Edit: Simpler question:
Why does this KMZ work:
https://drive.google.com/file/d/0B9n9YNu8ulsmb1ljQmJMdW5QQXc/edit?usp=sharing

But not this one:
https://drive.google.com/file/d/0B9n9YNu8ulsmLXYwNkpIQlY2STg/edit?usp=sharing

/Edit

Code:

import os
import zipfile
import arcpy

folder =r"C:\KMZTest"
table = r"C:\ProjectNumberSwitch_20140729.xlsx\C_Projects_Being_Switched$"
tblcolm = "Old_P6_NBR"
temploc = r"C:\KMZTest\Temp"

oldP6list = []
print "Finding project numbers"
cursor = arcpy.da.SearchCursor(table, tblcolm)
for row in cursor:
    oldP6list.append(row[0])
del row
del cursor

updatedlist = []
for root, dirs, files in os.walk(folder):
    for fil in files:
        if "Projects.kmz" in fil or "Points.kmz" in fil:
            print
            print fil
            print
            flpath = os.path.join (root, fil)
            zipfl = zipfile.ZipFile(flpath)
            foundoldP6 = False
            for content in zipfl.namelist():
                zipfl.extract(content, temploc)
                if content == "doc.kml":
                    doc = os.path.join(temploc, content)
                    os.rename(doc, doc + ".tmp")
                    destination = open(doc, "w")
                    source = open(doc + ".tmp", "r")
                    for line in source:
                        for oldP6 in oldP6list:
                            if oldP6 in line:
                                foundoldP6 = True
                                print "found:", oldP6
                                line = (line.replace(oldP6, oldP6.split("C")[0] + "VC" + oldP6.split("C")[1]))
                        destination.write(line)
                    source.close()
                    destination.close()
            zipfl.close()
            if foundoldP6 == True:
                print "Updating", fil
                updatedlist.append(fil)
                os.rename(flpath, flpath + ".tmp")
                NewZip = zipfile.ZipFile(flpath, "w")
                for newfl in os.listdir(temploc):
                    if not newfl.endswith(".tmp"):
                        addedfile = os.path.join (temploc, newfl)
                        NewZip.write(addedfile)
                os.remove(flpath + ".tmp")
            for oldfile in os.listdir(temploc):
                os.remove (os.path.join(temploc, oldfile))
print
print
print "Updated:"
for update in updatedlist:
    print update

print
print "Done"


texts: Creating image with specific latitude/longitude positions using GDAL?

texts:
I have an ASCII file with latitude, longitude, and data_val in the following format.

35-13.643782N, 080-57.190157W, 118.6
...

I have a GeoTiff image file, and I can easily view it.

I want to place a "pin" (can be a dot/flag/star or whatever is easiest) on the image at the specific latitude/longitude position found in the ASCII file.

Here is what I've managed to do so far:

My source image looks like this:

Driver: GTiff/GeoTIFF
Files: /tmp/Charlotte SEC 100.tif
Size is 16867, 12358
Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",38.66666666666666],
    PARAMETER["standard_parallel_2",33.33333333333334],
    PARAMETER["latitude_of_origin",34.11666666666667],
    PARAMETER["central_meridian",-78.75],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-365041.822331817995291,240536.419747152860509)
Pixel Size = (42.334586069440391,-42.334898968590878)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2016:06:24 12:46:45
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -365041.822,  240536.420) ( 82d48'55.43"W, 36d13' 4.92"N)
Lower Left  ( -365041.822, -282638.262) ( 82d35'10.11"W, 31d30'17.00"N)
Upper Right (  349015.641,  240536.420) ( 74d51'46.40"W, 36d13'26.16"N)
Lower Right (  349015.641, -282638.262) ( 75d 4'55.60"W, 31d30'36.99"N)
Center      (   -8013.091,  -21050.921) ( 78d50'12.11"W, 33d55'36.35"N)
Band 1 Block=16867x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 255,255,255,255
...

Here is what I've managed to cobble together in Python:

from osgeo import gdal, osr

src_filename = '/tmp/Charlotte SEC 100.tif'
dst_filename = '/tmp/foo.tiff'

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Specify raster location through geotransform array
# (upperleftx, scalex, skewx, upperlefty, skewy, scaley)
# Scale = size of one pixel in units of raster projection
# this example below assumes 100x100
gt = [-365041.822, 100, 0, 240536.420, 0, -100]

# Set location
dst_ds.SetGeoTransform(gt)

# Get raster projection
epsg = 4269            # http://spatialreference.org/ref/sr-org/lambert_conformal_conic_2sp/
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()

# Set projection
dst_ds.SetProjection(dest_wkt)

# Close files
dst_ds = None
src_ds = None

But, I can't quite figure out how to place a "red dot" at
35-13.643782N, 080-57.190157W

I'm having to learn some new details here (nomenclature about GIS).


texts: How to perform SQL query using .STContains in MS SQL Server

texts:
I have a table that stores polygons and a table that stores points with a value field. Both sets of data are stored using the Geometry type in MS SQL Server.

I have a simple requirement - I need to retrieve the records of points that are contained within the polygons and then sum the result based on a value field in the points table. So, the end result of the query is two columns:
- polygon names
- sum of the value field in the point table

DECLARE @area GEOMETRY  = (SELECT Shape FROM GREENFIELDAREAS WHERE AREANAME = 'Aerodrome')
SELECT NETDWELLINGS, Shape FROM ORA_NETDWELLINGCOMPLETIONS WHERE @area.STIntersects(Shape) = 1

The above query returns point geometries but only for 1 polygon (Aerodrome). What I need is to sum NETDWELLINGS and return that value for all the polygons.


texts:
I have a simple project — a base map and a gpx track. The basemap is from OpenLayers. I want to change the map the CRS to Arizona Central so I can use miles as units for QChainage.

When I change away from the default WGS84, my GPX track disappears.

I've tried to set the Layer CRS for the track. I've tried setting the CRS for the baselayer (the layer doesn't show), but then setting the CRS for the project from that layer (which at that time isn't showing), and then everything shows but the CRS has been reset to WGS84.

What am I doing wrong?

PS: I originally posted this at StackOverflow and saw another posting that suggested this was a better place to post.


texts: Calculating distance travelled along track in QGIS?

texts:
I'm new to GIS.

I have some tracks of people walking from point A to point B that I took using breadcrumbing with Garmin e-trex (set on the Auto track settings). I've loaded the tracks into QGIS and would like to calculate distance travelled along each track (cumulative distance is most important, but distance from trackpoint to trackpoint could also be useful). However, in the Attributes Table I can see only coordinates data, not distance data.

Is there any easy way to calculate distances travelled along tracks in QGIS?


texts:
I am trying to enable an enterprise database.
I right click > enable and then navigate to our organizations keycode file.
I get this error: Attribute Column Not found... (see screenshot)

This process has worked in the past but maybe IT has tweaked things or since 10.4 install this may have also changed things?

specs:

ArcGIS 10.4 (Advanced License)
SQL Server 2008
Microsoft SQL Server 2008 R2 (RTM) - 10.50.1600.1 (X64)


texts: Conflict with ArcGIS for Desktop and Erdas Emagine on same PC?

texts:
I already installed ArcGIS 10.1 for Desktop.

The problem is when I install th Erdas Imagine 8.6 that specify above.

So, could you help me in order to run both software products in one PC.


texts: ModelBuilder parameters not changing when running in batch?

texts: I have a simple model that selects layer by attribute of the feature class with the parameters of start year and end year set up as doubles and then use an expression in the select layer by attribute as: datesample_1 >= %Start Year% AND datesample_1 <= %End Year%.  It then uses copy feature to make an output feature class.  When I run the model individually it works fine, but when I run it as a batch even though I'm giving it different parameters, when it runs it uses the first parameters I entered over and over so my output features are all the same.


texts: Maximum angle between curved characters makes QGIS curved label invisible

texts:
When labeling a curved line in QGIS 2.2, and setting Placement to 'Curved', the behavior of 'maximum angle between curved characters' is not intuitive, and I haven't been able to find documentation on it.

The line I'm trying to label is pretty jagged (it's a hiking trail).  If I leave the 'maximum angle' settings both at 20 (default), no label is drawn.  If I crank the values up to 45, the label is drawn but it is very jagged just like the trail, and illegible.  If I crank the values down to 30, the label disappears again.

So, it seems like these values are being interpreted like so:
"If any generated character would be more than x degrees rotated from the previous character, then do not draw the label."

Whereas I was hoping the values would be used as a sort of filter / smoothing algorithm:
"Rotate the characters as needed such that each character's rotation is within x degrees of the adjacent characters' rotations."

Can anyone provide a more clear interpretation?


texts: Drawing line and buffer around it onto spheroid using PostGIS

texts:
I need to draw a line between two points and create a buffer around this line. The only difficulty is that I need to draw not a "cartesian" line but a line that is laying on the spheroid.

I tried to use ST_MakeLine but it has only geometry parameters so it draws only a straight line, not a curve (which is a straight line on spheroid).

Unfortunately ST_LineFromText also does the same.

Here is an example from http://postgis.net/workshops/postgis-intro/geography.html.

I want to draw the red curve. (Btw the length of this curve can be calculated by the ST_Distance(geography, geography))


texts: Extract soil and slope data using coordinates

texts:
I have no experience at all in GIS and I'm trying to extract soil type and slope data.
All I have is a CSV file with latitude and longitude coordinates.

Is it possible to use this to extract the soil type and slope data of each point?

I am using QGIS although I barely know anything.


texts: Access attribute names and add new attributes in processAlgorithm

texts:
How can I access the layer attribute names of the selected layer?
I have tried checking on Google but all the methods are mostly for the previous version and I am using py3 version.
att = source.fields()

but it does not return any results. I am unable to find a suitable solution for the problem.


texts: How to fill values in a field of a layer if the id columns of two layers match using pyqgis 3?

texts:
I have two layers, pipes and intersect layers. Both layers have an id, score column but intersect layer does not have all the ids present in pipes. So I am trying to check if the id present in pipes exists in intersect, if yes then copy intersect's score to pipe's score else leave it blank and continue. But the code neither gives an error nor does it make the changes in the attribute table of the pipe layer. Am I missing something?
The code is as follows:

pipe_feat = selectedLayer.getFeatures() #pipes layer
field4 = selectedLayer.fields().indexFromName(fid1) #id column is referenced
selectedLayer.startEditing()
featLayer = layer.getFeatures() #intersect layer
field3 = layer.fields().indexFromName('PIPEID')
for feat in pipe_feat:
    for f in featLayer:
        if feat[field4]==f[field3]:
             feat['SewerScore'] = f['Score']
             selectedLayer.updateFeature(feat)
        else:
             continue
selectedLayer.commitChanges()

I have also tried this method: But even this does not work.

  id_score={}
  for f in featLayer:
      id_score[f[field3]] = f['Score']

  for feat in pipe_feat:
      for key,value in id_score.items():
          if key == feat[field4] :
              feat['SewerScore'] = id_score.get[key]
              layer.updateFeature(f)
              selectedLayer.updateFeature(feat)
          else:
              continue


texts: latitude longitude into corresponding location in a image

texts:
I have an image , some data on this image, some data on the camera that took that image and a piece of code.
And i'm trying to find the location (Px,Py) in the image from (lat,lon).

my problem is that i don't know what is R0 , A1 , A2 from the camera data. how many values they have?

How do i calculate the radial dist. of A1 and A2?

And how do i calculate the distortion of R0?

(i got the camera data the way i posted it below.. messy)

I did mange to get x , y from collinearity equations using image rotation matrix, focallength... (after converting the data to ecef)
(lat,lon,alt) = (48.13564, 11.58147, 507)    // altitude from DEM

after calculation:    x = 0.0061  , y = 0.0523

Now i got this code that suppose to find the location (Px,Py) in the image using that (x, y) that i found:

    // x and y from collinearity equations using image rotation matrix, focallength …

 //orientation x in –row direction
    //orientation y in –col direction

    r = sqrt(x*x+y*y) ;
    dr = interior->radialdist_A1*(r*r-interior->distortion_R0*interior->distortion_R0)
                 + interior->radialdist_A2*(r*r*r*r-interior->distortion_R0*interior->distortion_R0*interior->distortion_R0*interior->distortion_R0) ;
    dx = dr*x + interior->xPrincipalPoint;
    dy = dr*y + interior->yPrincipalPoint;

    *x_out = (x + dx ) ;
    *y_out = (y + dy ) ;

    Row = -x/pixelsize + centerrow
    Col = -y/pixelsize + centercol

**

image data:

**

opk   3.99086  -2.62820 -68.23817

rpy  -1.28306   4.37035 110.28155

XYZ 692053.711 5334599.154  1612.905 (in utm zone 32)

LONLAT  11.58147  48.13564

Ground To Image Rotation Matrix

0.37035919   -0.92766421   -0.04767790

0.92775612    0.36688624    0.06828668

-0.04585474   -0.06952406    0.99652584

**

Camera data:

**

Col row  pixelsize                image center              R0            principal point xp, yp             A1                        A2                                                                                focallength

3744 5616  6.410000e-06  1872.000  2808.000   0.01400 6.793401e-07 2.887295e-05 -2.276840e+01 1.239733e+04 1.489102e+05 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 5.168871e-02


texts: Gridding or interpolating elevation data?

texts:
I wonder if someone can recommend a method for gridding elevation data in QGIS, GRASS or similar?

The elevation data I have is collected along 1m spaced lines at approximately 1 reading per second. So far I've tried inverse distance weighted, bicubic spline, and triangulation methods for the elevation data but all these produce a striping effect in the direction which the data is collected. I am not entirely sure that the error is due to the methods of gridding/interpolation I have been using. However, in tandem with the elevation data I am collecting a series of magnetic measurements (higher density), and these show no linear striping.


texts: How to force the urlUtils proxy perform only POST requests

texts:
How can I force the   urlUtils to use only POST requests.
I a using it as shown in the following block described also in the how to use proxy section, and I have noticed that some times (maybe when the request is less than 2000 characters ) it performs GET request, which in my case is not an accepted solution.

urlUtils.addProxyRule({
urlPrefix: "route.arcgis.com",
proxyUrl: "<url_to_proxy>"
});


texts:
Doing network analysis in QGIS. This is where I came from:

Start situation: + 600 points with fid and 1 relevant attribute. Also road network data.
First step: ran "network tools - service area (from layer)" with 3000 meters. Result: line features covering 3000 meters of the road network from each point.
Second step: the same thing but 5000 meters.

So in this point I have two line data layers, both with +600 features. One row in each table represents one service area and has the original fid + one relevant attribute.
Now I should subtract the 3000m service areas from the 5000 service areas in order to find those lines that are exactly between 3000m and 5000m from the point. How should this be done? I of course want to subtract the corresponding fids (eg. 5000m fid =145 subtracted with 3000m fid=145).


texts:
I'm try to load this shapefile with Basemap in order to plot it in Python:

http://geodata.gov.gr/en/dataset/periphereies-elladas/resource/7c80a2c1-93b7-4814-9fc4-245e775acaa6

But I get the following error:

ValueError: shapefile must have lat/lon vertices - it looks like this one has vertices
in map projection coordinates. You can convert the shapefile to geographic
coordinates using the shpproj utility from the shapelib tools

I tried using shpproj or ogr2ogr as another answer in this site suggested, but for the life of me I can not figure out what parameters I should pass to those programs, as I have zero knowledge about GIS in general. So if anyone here could help me with the conversion command that would be great.

I am on Ubuntu and the code I'm using is pretty simple I guess:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
map = Basemap()

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#ddaa66',lake_color='aqua')
map.drawcoastlines()

map.readshapefile('periphereies/periphereies', 'periphereies')

plt.show()


texts:
I'm looking for python documentation for specific operations in the Qgis python console. I know there is no suggestions simular to Arcpy but I would like a pure list of operation code lines as "layer.featureCount()"

Or is this taken from python standard librarys?


texts: QGIS Refactor Fields tool for altering attribute tables states: "This algorithm allows editing the structure of the attributes table of a vector layer. Fields can be modified in their type and name, using a fields mapping." Fields can also be added and deleted. However, output also alters the  length of all fields to 255. Is this expected behavior? Can original field lengths be retained?


texts:
Is it possible to show only every other contour line? And how do you do that?

(I'm new to QGIS Lisboa 1.8.0, but even newer to QGIS Dufour 2.0.1, but have used ArcGIS a bit a while ago, and also some ArcView 3.x earlier. But I'm now using QGIS Dufour 2.0.1. for this).


texts: Shortest distance from polygon centroids to a line in QGIS

texts:
To get the distance of each polygon centroid from a line, I have tried the Distance to nearest hub algorithm from the Processing Toolbox. The distances I got from the algorithm did not seem correct, upon checking manually. I later tried the v.distance GRASS command from within the Toolbox with the following settings:

It could be that I am not specifying the options correctly - for example, my "To" layer is a line and I have set the "Column where values specified by ' upload'..." to the "Descriptio" field. If this means that the calculated distances will be uploaded into the field then it could be the source of the issue. However, there is no option to create a new field to hold the distances (besides, a new layer is anyways created so I dont quite understand this setting). Nevertheless, running this command gives an error and QGIS stops working:

I dont know what the error is referring to. Is there some other plugin/algorithm that I can use? Alternatively, what are the correct settings for Distance to nearest hub and v.distance? I would like to achieve this entirely within QGIS, but am open to using R as well.

PS: I am running the latest version of QGIS (2.4 Chugiak) on Win 7 x64. I have already referred to the following questions on GIS SE:

Shortest distance from point to line (big dataset)

Measuring the distance between lines and points in QGIS

UPDATE: I could not get this to work from within QGIS despite trying several permutations of the GRASS v.distance settings. I finally used the gDistance function from the rgeos R package to obtain the result. However, while I have solved my problem for now, I would still like to know if one can do this is QGIS.


texts: Extracting LANDSAT time-series data by polygon using Google Earth Engine?

texts:
importing the landsat 5 TM 32-day NDVI composite data set.

I have made it as far as drawing a polygon of my ROI and extracting from a different data set, a global land cover map (ESA/GLOBCOVER_L4_200901_200912_V2_3) which was linked from another question posted on here:

Drawing polygon and extracting land cover inside using Google Earth Engine?

I tried to reconfigure it to work with my Landsat NDVI set and I was met with errors. I tried to extract just the NDVI band and try with that, but I keep getting an error:

  land.clip is not a function in , line #

var landsat: Imagecollection "Landsat 5 TM 32-Day NDVI Composite"
var geometry: Polygon, 4 vertices
var globcover: Image "GlobCover: GLobal Land Cover Map"

// Extract the landcover band
var land = globcover.select('landcover');
// Clip the image to the polygon geometry
var landcover_roi = land.clip(geometry);
// Add a map layer of the landcover clipped to the polygon.
Map.addLayer(landcover_roi);

reworking to accomodate the landsat:

// Extract the landcover band
    var land = landsat.select('NDVI');
    // Clip the image to the polygon geometry
    var landcover_roi = land.clip(geometry);
    // Add a map layer of the landcover clipped to the polygon.
    Map.addLayer(landcover_roi);

This is where I get the error,

land.clip is not a function in <global> line 4.

I am attempting to draw a polygon (square) around an area on the LS5NDVI data, export the "data cube", and correlate the time NDVI values at each pixel against various climate data to produce a single layer "heat map" of sorts to visually represent the correlation at each pixel.


texts: Python-arcpy : How make feature list from gdb and value list from field and use them in gaLayer function?

texts:
I'm not an expert in python, however I have to write a script to do geostatistical analyses with the GACreateLayer function from Arcgis.

I have several data (>2000) and I want to use list of features and list of field values  to perform cokriging.

I have written a script, the first part return list of features in my gdb and the list of values in the field that I have chosen (it works and return field values and feature names) and the second part do cokring (it works for one feature and its field)

However when I want to use both lists (inputDSet, end of script) in order to create a GA Layer it doesn't work.

File "C:\Program Files\ArcGIS\Desktop10.0\ArcPy\arcpy\ga.py", line 921, in GACreateGeostatisticalLayer     raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 045001: Input dataset(s) error. Table of inputs is not complete.
Failed to execute (GACreateGeostatisticalLayer)

I don't know how correctly use lists? May I transform them or something else?

Here, is my code:

import sys, string, os, arcgisscripting
# Create the Geoprocessor object
gp = arcgisscripting.create()

def listFcsInGDB():
    # list all Feature Classes in a geodatabase
    gp.workspace = "C:\\INTERP\\PJ_class.gdb"

    fcs = []
    for fds in gp.ListDatasets('','feature') + ['']:
        for fc in gp.ListFeatureClasses('','',fds):
            #yield os.path.join(fds, fc)
            fcs.append(os.path.join(fds, fc))
    return fcs

fcs = listFcsInGDB()
for fc in fcs:
    print fc

# List of value in fields
input_dataset =  fc
Atts = 'PJ_RacPJ'    #field with rainfall values

rows = gp.searchcursor(fc)
row = rows.next()

NewList = []
for fc in fcs:
    for row in gp.SearchCursor(fc):
        fcValue = row.getvalue(Atts)
        NewList.append(fcValue)
#print NewList

arcpy.CheckOutExtension("GeoStats")

# Load required toolboxes
arcpy.ImportToolbox("C:/Documents and Settings/ArcGIS/Toolbox.tbx")

# Local variables:
Cokriging_xml = "C:\\INTERP\\COK_.xml"
CK = "CK2"
ids1 = "C:\\INTERP\\data.gdb\\dem1000"

for fc in fcs:
    for fcValue in NewList:
        inputDset1 = fcValue          #Variables
        inputDset2 = ids1            #Covariable
        InputDset = "inputDset1;inputDset2"

# Create cokriging layer (code from model builder)

    tempEnvironment0 = gp.autoCommit
    arcpy.env.autoCommit = "1000"
    tempEnvironment1 = gp.spatialGrid1
    arcpy.env.spatialGrid1 = "0"

    arcpy.GACreateGeostatisticalLayer_ga(Cokriging_xml,InputDset,CK) ##Function

    arcpy.env.autoCommit = tempEnvironment0
    arcpy.env.spatialGrid1 = tempEnvironment1

I need help.

If someone has an idea?  It would be great.

Thanks


texts: g.getOutputs is not a function while Creating a Chart from a Series

texts:
I created a Time Series of NDVI from Sentinel 2 Data, now I want to Chart it, using the following code. I get the error

  g.getOutputs is not a function

I never used g.getOutputs, and from what I read, it seems to be an error on the GEE side.

 var sentinel2Collection = ee.ImageCollection('COPERNICUS/S2')
 .filterDate('2000-01-01','2017-12-31')
 .filterBounds(aoi);

 print(sentinel2Collection);

Map.centerObject(aoi, 10)

var imageS2 = sentinel2Collection.select(['B8', 'B4']).mosaic();

var addNDVI = function(im){
    var ndvi = im.select(['B8','B4']).normalizedDifference();
    return ndvi;
};

var ndviCollectionS2 = sentinel2Collection.filterBounds(aoi).map(addNDVI);
print(ndviCollectionS2.first());
Map.addLayer(ndviCollectionS2,  { min: -0.8, max: 0.8,'palette': 'FF0000,FFFF00,00FF00'},' NDVISentinel2');
print(ndviCollectionS2)
ui.Chart.image.series(ndviCollectionS2, /*'Sentinel2NDVi',*/ aoi, 'median');

Code to GEE


texts: Tiff file output using print composer has blocks of image scrambled (out of place)

texts:
I have used 2.18 to print maps to file (tiff) using print commposer.  The base layer imported into QGIS consists of NOAA raster tiles. (My problem arises whether or not I have added other layers.)  The composed maps seen in the composer look fine.  However, the tiff file(s) created have problems.  In particular, the tiff file output shows "sub blocks" of the image in composer moved from their original locations (e.g., scrambled).
My sense is that the problem arises from my use (misuse?) of raster tiles.

Any ideas regarding the nature of the problem and solutions/fixes?


texts: How can I select features with both multiple attributes and NULL values

texts:
I'm using ArcMap and am trying to filter a selection from my attribute table that contains both NULL and regular values but it is not working.

The query should only select John values with NULL, Addition, and Construction. Below I did my best but the NULL are not helping.  What is the solution?

Select * FROM Table_name WHERE:
Usr_name ='John' AND Project_name IS NULL OR Project_name IN ('Addition', 'Construction')

Input

Object_ID   Object_type Usr_name    Project_name
 1            Triangle    John         Addition
 2             Square     Peter        Addition
 3            Polygon     Mary         Addition
 4             Circle     John       Substracion
 5            Polygon     Mary          <Null>
 6             Square     Mary          <Null>
 7             Square     John          <Null>
 8            Polygon     John         Obstacle
 9              Hex       John       Construction

The minimal desired output result

Object_ID   Object_type Usr_name    Project_name
  1           Triangle    John        Addition
  7            Square     John         <Null>
  9             Hex       John      Construction


texts: Mosaicking a Image Collection by Date (day) in Google Earth Engine

texts:
I have an image collection of Sentinel-1 data over 3 years filtered by a region of interest that spans 3-4 tiles, and want to know the best way to mosaic all the images (tiles) that were taken on the same day into one image. Currently filtering by date yields separate images for each separate tile, and I would like to merge all images (tiles) from the same day into one image.

I know how to do this for one specified date with a command like this:

ee.Image(collection.filterDate('2014-10-01','2014-10-02').mosaic())

But I would like to do this iteratively for the whole collection in my ROI over a broad date range and obtain a new image collection of mosaicked images by days in which there is some data.

The collection is created by the following code:

var poly = /* color: #d63000 */ee.Geometry.Polygon(
    [[[-76.0803, 10.8656],
      [-76.0913, 7.7436],
      [-73.1909, 7.7545],
      [-73.3776, 9.4273],
      [-75.2124, 10.9304]]])

var start = ee.Date('2014-10-01');
var finish = ee.Date('2018-03-31');

var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterDate(start, finish)
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filterMetadata('resolution_meters', 'equals', 10)
.filterBounds(poly);

For example, the first three images in this collection are all tiles from the same day (2014-12-02), and I would like to mosaic them together into one 3-tile image for day 2014-12-02, and then do this for the entire  collection of over 500 tiles/images to get a collection of day-unique tile mosaics.


texts: Printing maps at a specific scale using the print composer (1"=2000')

texts:
I have created my map for land planning and set up my sheet to be exported in the print composer, but I'm confused when it comes to printing at a certain scale. I see within the item properties > main properties that there is a tab for "scale" but I have no idea what unit of measurement it is displaying at.

I have attached images of my coordinate system and general settings for reference, and in the case of this project I am looking to create a basemap at a scale of 1"=2000'.


texts: Using thematic layer to restrict segmentation and classification in eCognition

texts:
This question has been asked before as Using thematic layer for segmentation in eCognition?.
However, I have not been able to resolve the issue using the answer provided.
Can someone clarify the correct procedure?
I would like to restrict all segmentation and classification processes to take place within areas (polygons) covered by my thematic layer rather than the whole scene.
I have successfully carried out a chessboard segmentation using my thematic layer so the polygons from my thematic layer are now image objects. I then used the 'assign class' algorithm to assign those image objects to one class. I then attempted to carry out an additional segmentation at the image object level using the class filter to select my new class. However, the segmentation was still carried out on the entire scene.
The only instance where my steps have differed to those recommended at Using thematic layer for segmentation in eCognition? is when I was unable to see how to assign a class based on the "number of overlaps". I assigned my image objects to a single class without using "number of overlaps".


texts:
I'm trying to create a multipolygon dataset where each row is identifiable.

Example below

library(sf)
library(raster)

# example raster data
r <- brick(nl = 2, nrows = 5, ncols = 10)
set.seed(54); r[] <- rpois(ncell(r)*nlayers(r), lambda = 3)
names(r) <- c("layer_1", "layer_2")
r

for (i in 1:nlayers(r)) {
  s <- r[[i]]
  # set values below x to NA, all else to 1
  s[s[] < quantile(r[[i]], 0.9)] <- NA
  s[!is.na(s[])] <- 1
  # convert from raster to multipolygon
  s <- st_union(as(rasterToPolygons(s, na.rm = TRUE), "sf"))
  # send to workspace
  ## ideally here I could bind to an existing empty dataset
  ## after naming to rows
  assign(x = paste0("poly", i), value = s)
}
# no way to tell which one is poly1 or poly2?
poly3 <- c(
  st_cast(poly1, "GEOMETRY"),
  st_cast(poly2, "GEOMETRY"))
poly3
names(poly3) <- names(r)
names(poly3)
write_sf(poly3, "./test_poly_output.shp") ## works
rownames(poly3) <- names(r) ## error
poly3$name <- names(r)
poly3
write_sf(poly3, "./test_poly_output_2.shp") ## error

I want to be able to open the dataset in QGIS or Arc, and see which rows correspond to the appropriate layer. Opening the test_poly_output.shp shows the following which isn't particularly helpful, even though in this example FID_1 == 0 obviously corresponds with layer 1

Ideally I would be able to bind to an existing empty sfc class. The code below works for a MULTILINESTRING, but I cannot get it to work for MULTIPOLYGON or GEOMETRY types.

## Example where LINESTRINGS can be appended to blank sfc
rConts <- st_sf(st_sfc(crs = 4326))

for (i in 1:nlayers(r)) {
  contsR <- r[[i]]
  contsR <- rasterToContour(r[[i]], levels = quantile(r[[i]], 0.9))
  names(contsR) <- "Q90"
  contsR <- as(contsR, "sf")
  rownames(contsR) <- gsub("X_", "", names(r[[i]]))
  contsR$Year <- as.integer(gsub("layer_", "", names(r[[i]])))
  rConts <- rbind(rConts, contsR)
}
rConts
write_sf(rConts, "./test_contourLines_output.shp")


texts: GRASS GIS Addon Installation Windows 10 Problem

texts:
I know there are some very similar questions regarding this topic:
GRASS GIS addon "extension not found"; could the problem be with the GRASS Addons server?

but so far nobody has found a proper solution it seems. I am trying to install r.stream.distance (and other) addons but not with Grass 7.6 or 7.8 i get any success. I tried to run a different g.extensions script (as suggested in the question mentioned above) and the addon seems to install but then there are alwas problems in running the addon, specifically:


texts:
I have a shapefile that has 2 attributes named "X" and "Y".  I want to change this to "Lon" and "Lat".

My (limited) experience tells me that I should read in the original shapefile, preserve the original driver, crs, make a new schema with "Lon" and "Lat" defined, and then write a new/new named shapefile.

Is there a quicker/easier way with fewer steps?

import fiona

with fiona.open(filein, 'r') as source:
    source_driver=source.driver
    source_crs=source.crs
    source_schema=source.schema
    print(source_schema)

new_schema= {'geometry': 'Point',
         'properties': { 'Lon': 'float:16.6',
                         'Lat': 'float:16.6',
                         'SP': 'float:16.6',
                         'X_meter': 'float:19.11',
                         'Y_meter': 'float:19.11',
                         'Depth': 'float:19.8',
                         'Line_Name': 'str:50',
                         'SurveyID': 'str:50'}}

with fiona.open(fileout, 'w',
        driver=source_driver,
        crs=source_crs,
        schema=new_schema) as c:
    print(len(c))

c.closed


texts: How to calculate for every basin/sub-basin, drainage area at every elevation contour along stream?

texts:
With DEM as a source, I created all necessary files like flow direction, accumulation, slope, streams, watersheds, pour point etc... In my watershed (basin), I have many sub-watersheds (sub-basin) with streams and outlets at various elevations.
I need to do analysis on the following two basis.

I need to calculate for every basin/sub-basin, the drainage area at every elevation contour along the stream. Please ref figure for better understanding my view.
I the figure there are six elevation contours (2 black top-bottom, 4 red in between) in the basin/sub-basin. In the figure, outlet is at the bottom elevation. Now I want to change outlet to all above elevations along the stream. And for each elevation point outlet along the stream, I need to calculate drainage area and volume.

I am trying various options with ArcGIS, Global Mapper, AutoCAD, Surfer, MicroStation and others...

Please give me some idea in doing this watershed analysis with less effort with any of the above said software/tools/extensions/toolbox with better accuracy.


texts: I want to create a fly-through in GRASS GIS 3D. The problem: throughout the "flight" the rendering is at low visual resolution, while it is set a high resolution. Is there a way to render a flight sequence at high resolution (1)?


texts: Single-ringed GeoJSON polygon in Mongodb

texts:
I've loaded the Census' zip code data into a mongo database.

I can query specific zip codes like this:

db.zips.findOne({"properties.ZCTA5CE10":"80205"})
db.zips.findOne({"properties.ZCTA5CE10":{$in:["80205", "80202"]}})

For the second mongo find, this is what those zip codes would look like on a map:

Now I want to query specific zip codes and "merge" any that are tangent creating single polygons.  For example, in the second example I want to draw those two polygons as a single polygon.

Is there a way to do that with mongodb?

In postgis I think I could accomplish this with ST_Collect or ST_Union (related stack overflow question).


texts: Is there a possibility to simply display the next element on a layer?

texts: I have to check all the elements on a layer and I'm looking for an easy way to go through all of them. One option would be to open the attribute table and just always click on the next line and then "move to selection". I would prefer an option where I just click on an arrow or use a shortcut and QGIS just jumps to the next element, so that I don't have to be aware of not missing a line and just needing half the amount of clicks. Is there anything like that in QGIS?


texts: How to get the layer name of the WMS in the legend of the user interface

texts:
I have a problem with displaying the name of the layer in the legend of the user interface. I prepare the mapfile and I access with the WMS in my local test interface. When I add it to the user interface, the WMS appears but, the name of the layer does not appear in the caption. Here's my mapfile code
LAYER
    INCLUDE "../shared/ec/ec_getFeatureInfo.map"

    NAME "WMS"
    GROUP "Bâtiments"

    METADATA
        INCLUDE                 "../shared/ec/ec_meta_layer.map"
        "wms_group_title"       "Bâtiments"
        "wms_title"             "Extraction automatique d'empreintes de bâtiments"
        "wms_name"              "automatic_extraction_building_footprint_fr"
        "wms_server_version"    "1.3.0"
        "wms_format"            "image/png"
        "wms_dataurl_format"    "igo"
        "wms_dataurl_href"      "https://ouvert.canada.ca/data/fr/dataset/7a5cda52-c7df-427f-9ced-26f19a8a64d6"
        "wms_dataurl_format"    "text/html"
    END

    #MAXSCALE 10000000

    CONNECTION "http://cartes.geogratis.gc.ca/wms/automatic_extraction_building_footprint_fr?"
    CONNECTIONTYPE WMS
    TYPE RASTER

    PROJECTION
        "init=epsg:4269"
    END
END

here the display of the WMS to test it on my local server, in the legend the name of the layer appears

However in the user interface, it does not appear.
There would be a lack of functionality in the mapfile for the name of the layer to appear in the legend?


texts:
I have a large shapefile (1 gb) and need to rasterise it. I have tried the following already.

1.) Import it into GRASS using v.in.ogr this failed with the error message: ERROR: G_realloc: unable to allocate 498240036 bytes at break_polygons.c:188

2.) My second idea was to use PostGIS. Import the shapefile, resample it at the x,y location of the grid and then export these points and create a grid from xyz. I successfully imported the shapefiles (polygons and points) but intersecting 1 million polygons with 300k points seems to be very slow. I used the following PostGIS satement, maybe there is room for improvement.

select polygons.land_id,grid.geom from grid,polygons where grid.geom && polygons.geom and within(grid.geom,polygons.geom)

3.) I did also try to use simplify() in PostGIS. But I lost to many small polygons (i.e. some areas that were covered only with small polygons became null).

Any ideas would be greatly appreciated.


texts: QGIS 3.12 - How to divide Polygon even areas

texts:
I'm trying to divide a 5ha paddock into .5ha "cells".

In QGIS terms I'm looking to subdivide an irregular polygon into sub-polygons which will have the same area inside and as near as possible straight dividing fences. Ideally this would be two steps:

Divide 5ha paddock into 2 2.5ha thin paddocks with a vertical line (red line below)
Divide each 2.5ha paddock into five 0.5ha paddocks with horizontal lines (blue lines, should be straight.

Ideal after result:

I have seen the plugins Polygon splitter, divide polygon, and Split feature on steriods, but I cannot get any of them to work. At the end of this exercise I want a GPS route for the new fencelines


texts: Working with topology errors sample code error

texts:
I am following the sample code from ESRI for checking topology errors using vb.net. But in my workspace an error happened. My code is listed below:

   ' Step 1 - access a topology
    Dim topoUiD As UID = New UIDClass
    topoUiD.Value = "esriEditorExt.TopologyExtension"
    Dim m_application As IApplication = My.ArcMap.Application
    Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD), ITopologyExtension)  'original code from esri
    Dim topology As ITopology
    topology = CType(topologyExt.CurrentTopology, ITopology)

    ' Step 2 - Finding topology errors
    Dim geoDS As IGeoDataset = CType(topology, IGeoDataset)
    Dim errorContainer As IErrorFeatureContainer = CType(topology, IErrorFeatureContainer)
    Dim eErrorFeat As IEnumTopologyErrorFeature
    eErrorFeat = errorContainer.ErrorFeaturesByRuleType(geoDS.SpatialReference, esriTopologyRuleType.esriTRTLineCoveredByLineClass, Nothing, True, False)
    Dim topoError As ITopologyErrorFeature

    topoError = eErrorFeat.Next

    ' Step 3 - selecting topology errors
    topologyExt.ClearActiveErrors(esriTEEventHint.esriTENone)
    topologyExt.DelayEvents(True)

    While Not ((topoError = eErrorFeat.Next) Is Nothing)
        topologyExt.AddActiveError(topoError, esriTEEventHint.esriTENone)
    End While

    topologyExt.DelayEvents(False)

In the while sentence:  While Not ((topoError = eErrorFeat.Next) Is Nothing)

An error indicates:

Operator '=' is not defined for types 'ESRI.ArcGIS.Geodatabase.ITopologyErrorFeature' and 'ESRI.Geodatabase.ITopologyErrorFeature'. Use 'Is' operator to compare two reference types.

When I change "=" to "Is", another error indicates:

'Is' operator does not accept operands of Type 'Boolean'. Operands must be reference or nullable types.

So I don't know how to define the equals operator of 2 'ESRI.ArcGIS.Geodatabase.ITopologyErrorFeature'. Thanks in advance.


texts:
I'm trying to get a map of Manhattan that is aligned vertically. I have a GeoJSON file with polygons in EPSG:4326 in fname that I've been transforming to EPSG:6539

gdf = gpd.read_file(fname).to_crs("EPSG:6539")

Everything works as intended, but I would really like something like this example:

Is there a named transform that accomplishes this? If not, how to correctly rotate the whole CRS in GeoPandas? If I do something like gdf.rotate(28), I get each polygon rotated individually but not the whole map rotated, like so:


texts: Unwanted changes in layout of data-driven pages during export to JPG by ArcPy?

texts:
I'm using a python script to export a set of maps created with ArcMap's Data Driven Pages (DDP). Everything was working well until I added a new layer to the map - it automatically added to the legend, which I didn't want, so I removed it. Ever since then, however, the legend has been re-sizing at random during export.

It's not an overwrite issue (I'm removing all images from the target folder before running the script) and I've tried deleting and re-creating the legend, which has caused fewer of the maps to have re-sized legends, but there are still some. I'm not sure what to try next, or what caused this in the first place.

Here's my script (modified from Changing output name when exporting data driven pages to JPG?):

# Modify the directory to match your file locations
import_path_Singles = r"C:\Users\Keegan\Desktop\14Feb2017_RVCAmapping\Arc_v10_3\RMP_Ag_SinglePlates.mxd"

# The script will then generate output JPG files at a specified out-path:
export_path = r"C:\Users\Keegan\Desktop\14Feb2017_RVCAmapping\Arc_v10_3\OutputJPGs"

# Then specify a "sort" field
field_name = "ASSESSMENT_ROLL_NUMBER" # Name of field used to sort DDP
#field_name = "Assessment Roll Number"
print field_name

import arcpy
mxd = arcpy.mapping.MapDocument(import_path_Singles)
print mxd.dataDrivenPages.getPageIDFromName(field_name)
for i in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = i
    row = mxd.dataDrivenPages.pageRow
    arcpy.mapping.ExportToJPEG(mxd, export_path + "\\" + row.getValue(field_name) + "SinglePlate.jpg")
del mxd

The results can be seen in the image below - at left, good legend, at right, bad legend:

In Layout view, the legend looks like the one at left, and of the 21 maps I'm generating, some look like that one, and some look like the one at right.

Has anyone come across this before? Any idea what's going on here, and how to fix it?


texts: Geopandas errors when using ".within()" in an if statement

texts:
I have three dataframes I am working with and want to do the following logic:

if dataframe1 in not within dataframe2 but is within dataframe3 then do something

first attempt was:

if df1['geometry'].within(df2['geometry']) == False & df1['geometry'].within(df3['geometry']) == True:
    continue with code

and I got the error:

The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all() using panda python

so I tried:

if df1['geometry'].within(df2['geometry']).any() == False & df1['geometry'].within(df3['geometry']).any() == True:
    continue with code

and I get the error:

TypeError: 'GeoDataFrame' objects are mutable, thus they cannot be hashed

I'm not sure what to do to fix this

update:

tried

if df1['geometry'].any().within(df2['geometry']).any() == False & df1['geometry'].any().within(df3['geometry']).any() == True:
    continue with code

and I get the error:

AttributeError: 'GeoSeries' object has no attribute '_geom'


texts:
I am working on a publication and would like to make my schematic visually pleasing. The amount of detail is too high and the image would be better If I could simplify and smooth my line and polygon vectors. Ideally it would minic smooth Bezier like curves. (Curves are not yet supported in QGIS so searching an alternative)

I've tried many many settings in QGIS, primarily the v.generalize tool but I either get errors or the results still has sharp corners.

A close up reveals the ugly sharp corners of the vectors.

And I would like to get something like this:

Is there a simple way to simplify the vectors into smooth lines? Alternativily, which settings in GRASS v.generalize should I use?

Related question: Is there a way to render lines with a smooth (bezier) curves in QGIS?

It seems like v.generalize with Chaiken algorithm could do the trick. Searching for the correct settings. (Trial and error)

Shapefiles can be found here


texts: Showing ESRI Historical Topographic maps in Openlayers

texts:
I'm trying to show a Dutch historical map (see www.topotijdreis.nl) in Openlayers. However when I load the map it tries to load the tiles but the all give a 404 error. I'm guessing it has something to do with the Projection but I can't figure it out. Can someone point me in the right direction?

map = new ol.Map({
    target: 'map',
    layers: [
    new ol.layer.Tile({
        source: new ol.source.OSM()
    })
,
new ol.layer.Tile({
                   title: "1972",
                    type: "base",
                   source: new ol.source.XYZ({
                       url: 'http://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_1972/MapServer/tile/{z}/{x}/{-y}'
                   })
                })

],
    view: new ol.View({
        center: [592337.6034404513, 6824935.041410904],
        zoom: 7
    })
});

https://jsfiddle.net/r9xqf86e/18/


texts: Calculating Slope for the Whole World in PostGIS

texts:
The Challenge

I have a PostGIS Database which contains all of SRTM v4.1. I need to calculate slope for the entire globe as part of an error modeling process. The DEM is tiled at 100x100.

Whenever I try to calculate slope on the tiled data, I get edge effects. I've written about it and how to avoid it for small areas, but  I can't just ST_Union the globe and expect it all to be fine. Even for just the state of Colorado I get a memory allocation error.

The Plan

For each tile in the table, calculate the bounding box.
Expand that bounding box by 1 pixel on each side.
Using ST_Intersects, select all tiles that intersect the bounding box.
Union that 9-tile set, and calculate the slope.
Clip the resulting much larger raster with the original 100x100 bounding box
Insert my 100x100 slope tile into a new table.

Attempted Implementations

First Try

Based on tilt's answer, I've tried this query:

INSERT INTO tmp_slope(rast)
SELECT ST_Clip(ST_Slope(ST_Union(b.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(a.rast))
FROM dem_foco_100 a
INNER JOIN dem_foco_100 b
ON (ST_Intersects(a.rast, b.rast))
GROUP BY a.rast;

Unfortunately, as you can see from the links in my comments on his answer, this is still producing edge effects for me. It seems that this is the case because ST_Intersects doesn't return all neighboring tiles as expected.

Expanding the Bounding Box

Instead, I've created a query based on expanding the bounding box of the current tile with the JOIN route proposed by tilt to come up with something that works for my subset.

INSERT INTO tmp_slope(rast)
SELECT ST_Clip(ST_Slope(ST_Union(b.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(a.rast)) s
FROM (SELECT * FROM dem_foco_100) a
INNER JOIN dem_foco_100 b
ON ST_Intersects(
  ST_Envelope(b.rast),  -- All raster tiles
  ST_MakeEnvelope(      -- Expanded bouding box
    -- WEST
    (ST_UpperLeftX(a.rast) - ST_PixelWidth(a.rast)),
    -- SOUTH
    (ST_UpperLeftY(a.rast) - ((ST_Height(a.rast) + 1) * ST_PixelHeight(a.rast))),
    -- EAST
    (ST_UpperLeftX(a.rast) + ((ST_Width(a.rast) + 1) * ST_PixelWidth(a.rast))),
    -- NORTH
    (ST_UpperLeftY(a.rast) + ST_PixelHeight(a.rast)),
    -- Use the OG SRID
    ST_SRID(a.rast)
  )
)
GROUP BY a.rast;

The Cursor Approach

I've also tried a cursor using the same bounding box logic:

CREATE OR REPLACE FUNCTION slope_cx()
RETURNS void AS
$BODY$
DECLARE
    data CURSOR FOR select * from gdem_100;
    data_row record;
BEGIN
    OPEN data;

    LOOP
    FETCH data INTO data_row;
        EXIT WHEN NOT FOUND;

        INSERT INTO cur_slope(rast)
        SELECT ST_Clip(ST_Slope(ST_Union(a.rast), 1, '32BF', 'DEGREES', 111120), ST_Envelope(data_row.rast))
        FROM gdem_100 a
        WHERE ST_Intersects(
          ST_Envelope(a.rast),  -- All raster tiles
          ST_MakeEnvelope(      -- Expanded bouding box
            -- WEST
            (ST_UpperLeftX(data_row.rast) - ST_PixelWidth(data_row.rast)),
            -- SOUTH
            (ST_UpperLeftY(data_row.rast) - ((ST_Height(data_row.rast) + 1) * ST_PixelHeight(data_row.rast))),
            -- EAST
            (ST_UpperLeftX(data_row.rast) + ((ST_Width(data_row.rast) + 1) * ST_PixelWidth(data_row.rast))),
            -- NORTH
            (ST_UpperLeftY(data_row.rast) + ST_PixelHeight(data_row.rast)),
            -- Use the OG SRID
            ST_SRID(data_row.rast)
          )
        );

    END LOOP;
    CLOSE data;
END;
$BODY$
LANGUAGE plpgsql VOLATILE;

Both of these queries do what I want for small areas (< 600 tiles). Unfortunately, both fail on my full test area (~256000 tiles). Neither actually errored out, but I cancelled them after about 40 hours.

I need to cover the full area in an hour or less of run time for this to be feasible to run on the full GDEM table, as the land area of my full test dataset is ~1/90th the area of North America alone.

How can I write a performant query or function to calculate edge-less slope for the entire world?


texts: Color in QGIS based on record of events in CSV file?

texts:
I have a shapefile with the following attributes:

ID    Region
1     Region1
2     Region2
3     Region3

Also I have the following record of events in a CSV file

EventID     Year     Month     RegionID     Quantity
001         2015     12        1            6
002         2015     12        2            7
003         2015     12        3            3
004         2015     11        1            4
005         2015     11        3            3
006         2015     10        2            6
007         2014     12        3            3
008         2014     11        2            2
009         2014     10        1            7
...         ...      ...       ...          ...

I would like to color my map acording to a filtered record. For instance, say I want to display all the events of November 2015. Then I would have:

ID    Region    Quantity
1     Region1   4
2     Region2   0
3     Region3   3

Or say I would like to display all the events of 2015. Then I would have the sum of the events:

ID    Region    Quantity
1     Region1   10
2     Region2   13
3     Region3   6

Is there a way to solve this problem dynamically, without rewriting the attribute table of the shp every time?

Consider that the CSV file is updated frequently.


texts: Can I offset the map frame in the composer?

texts:
I have a map frame in a layout in which the map fills the frame. If I check the frame property checkbox, there is no whitespace between the contents and the frame. Is there a way to offset the frame? I know I can draw a rectangle around the frame and center the objects, but it's hard to get the aspect ratios identical.
qgis 2.01 64 bit
Windows 7


texts: Is it possible to bring PostGIS database Query results into QGIS workspace other than using DB Manager?

texts:
I would like to use the QGIS modeller (or similar) to import SQL query results form my Database directly into my workspace.

I am asking as part of my job is to repetitively add vector and DTM layers from our Server and clip them to a smaller size based on the Site boundary provided by the client. I have now managed to at least add the vector layers correctly, Haven't had much luck with the Raster layers yet TBH, and buffer this Boundary by 1km and insert the Layer using DB manager.

The SQL I have used is as follows-

SELECT vectors_index__21.gid, vectors_index__21.fid,vectors_index__21.geom, st_contains(ST_Buffer(public.site_boundary.geom,1000)public.vectors_index__21.geom)

FROM public.site_boundary, public.vectors_index__21

WHERE st_contains(ST_Buffer(public.site_boundary.geom,1000), public.vectors_index__21.geom) IS TRUE

GROUP BY vectors_index__21.gid,vectors_index__21.fid,vectors_index__21.geom, site_boundary.geom

Can I add this code to a QGIS model or script that will do this task multiple times, ie. with different vector and raster layers so that the workspace basically sets itself up each time once the new site boundary has been added to my database and SQL syntax updated?

I have looked eveywhere for this type of functionality but so far my only option seems to be the DB manager to import my results to a workspace. Has anyone done this before or have any suggestions?Everyone loves a good bit of automation in their lives.


texts: To what extent is SQL's language accessible via ArcGIS tools?

texts:
That question just sounds so redundant...
At any rate, I have been asked to explore the uses of SQL in ArcGIS. We use an SDE DB for all server-side storage and DB managing/large-scale editing. However, I am attempting to delve deeper into just how integrated SQL is in ArcGIS's tools. There is a note in the SQL Query Reference Guide of ArcGIS 10's help that clearly states SQL syntax does not work when using the Field Calculator tool. I find this very disappointing as it would seem very useful to be able to pull and mangle information directly from the database itself when calculating field values.

What are the largest limitations and pitfalls one might fall into when attempting to make use of SQL in ArcGIS? Which tools see best functionality when integrated with SQL?

If an example is required to make this a clearer question, I am attempting to track the number of points of interest that are of a certain type in a given region. Each point has an attribute which describes its type and region. My current solution involves heavy use of cursors for counting and data retrieval, but a superior pointed out that a simple SQL query could do all of this and more. Unfortunately, it seems the functionality is not present for this.

Or is it?


texts:
I would like to use the "create categorized renderer from style" Plugin in QGIS 3.14.15 but I didn't have any .XML-file for this.
I tried to create one for a vector layer but in the options to export my style I didn't find an XML extension. I exported the style to XLSX and converted it to XML. Then, I tried to import this with "create categorized renderer from style" and I set the input layer and the column for matching the data from the XML and set the style database (.XML-file) and ran it. However, I have been getting an error message:

Incorrect root tag in style: workbook.

How could I export an XML file to be used to "create categorized renderer from style" Plugin? I have to use this plugin for an assignment.


texts:
I have a set of time series daily raster files, I want to have a temporal combination on daily bases. As an example I want to combine day1 with day2 and day2 with day3 and so on.
The file name are exactly same and just date is varying in the file in the sequence of (2010001,2010002,2010003....).
I want to edit [if item1.endswith(item2[32:51]):] to match file with its proceedings file.

import  arcpy
import os
ws = arcpy.env.workspace = r'D:\BRB Snow Cover\Sample\Terra'
Terra_Rasters = arcpy.ListRasters("MOD*")
matches = 0
for item1 in Terra_Rasters:
    for item2 in Terra_Rasters:
        print item1
        print item2
        if item1.endswith(item2[32:37]):
        print "Match"
        output = "D:\\BRB Snow Cover\\Sample\\output""\\" + str(item1) + '.tif'
            arcpy.gp.Con_sa(item1,item2,output,item1,"\"Value\" =8")
        else:
            print "No"
            matches = matches+1


texts: Does MapInfo License Server run on Windows Server 2012 R2?

texts:
The documentation for MapInfo License Server Utility (v4.5) advises that it has been tested on Windows Server 2012.

I need to know whether the LSU runs on Windows Server 2012 R2


texts: Why is the arcgis.geometry.buffer method failing on supposed empty geometries?

texts:
I am attempting to buffer some output after querying a feature service using the arcgis API for Python:

from arcgis.gis import GIS
from arcgis.features import FeatureLayerCollection
from arcgis.geometry import Geometry, buffer

# obtain the LOMR feature
nfhl_url = "https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer"
nfhl = FeatureLayerCollection(nfhl_url)
lomr = nfhl.layers[1]

# define the desired spatial ref for the project
sr = 2876 # NAD83(HARN) / Colorado North (ftUS)

# start an anonymous GIS session
anon_gis = GIS()

# query the service
boulder_lomrs = lomr.query(where="DFIRM_ID = '08013C'", out_sr=sr)

# create a list of geometries to buffer
lomr_geoms = [l.geometry for l in boulder_lomrs.features]

# buffer the polygons
lomr_buffer = buffer(geometries=lomr_geoms,
                     in_sr=sr,
                     distances=[1],
                     unit='ftUS',
                     out_sr=sr,
                     buffer_sr=sr,
                     union_results=True,
                     gis=anon_gis)

The error I receive is:

Exception: Unable to complete operation.
The operation was attempted on an empty geometry.
(Error Code: 400)

However, I've checked that all geometries in lomr_geoms are valid and not-empty. This leads me to believe I'm doing something wrong with the spatial references and/or units.

Why am I getting this error? The ESRI documentation on the buffer method is not very specific on syntax, and they wrongly claim that unit will be derived from the in_sr parameter (I tried running the buffer without the unit parameter and it failed because it's required).

EDIT 1

I've tried redefining how the list of geometries is constructed:

# create a list of geometries to buffer
lomr_geoms = [Geometry(l.geometry) for l in boulder_lomrs.features]

# buffer the polygons
lomr_buffer = buffer(geometries=lomr_geoms,
                     in_sr=sr,
                     distances=[1],
                     unit='ftUS',
                     out_sr=sr,
                     buffer_sr=sr,
                     union_results=True,
                     gis=anon_gis)

I still get the same error code, but for a different reason:

Exception: Unable to complete operation.
(Error Code: 400)


texts:
I tried to manually set different scales for different features in atlas generation (QGIS 2.18.16). Thus, I had to switch off map - "controlled via atlas".

The resulting scales were fine, however the shapefiles were not centered anymore, but partly off the map. Any idea how to fix this?

As a workaround I now use the project-wide scales, but they are far from ideal. In the preview I could edit the scale, but it is apparently not saved.


texts:
The db structure is more or less like this and can't be modified. I can only add some indexes but can't set for example foreign keys:

    CREATE TABLE adresspoint
(
  id character varying(3999) NOT NULL,
  "position" geometry(PointZ,25832)
  CONSTRAINT pk_adresspoint PRIMARY KEY (id)
)

CREATE INDEX adresspoint_idx_position
  ON adresspoint
  USING gist
  ("position");

CREATE TABLE house
(
  id character varying(3999) NOT NULL,
  adresspointid character varying(3999),
  CONSTRAINT pk_house PRIMARY KEY (id)
)

CREATE INDEX house_idx_adresspointid
  ON house
  USING btree
  (adresspointid COLLATE pg_catalog."default");

I have query like that:

SELECT DISTINCT house.id
FROM house
   LEFT JOIN adresspoint ON adresspoint.id = house.adresspointid
             AND st_contains(st_makeenvelope('562827','6326412','563365','6326590',25832),adresspoint.position)
OFFSET 0 LIMIT 100

and it runs very quickly (100 ms) but when I change it to:

SELECT DISTINCT house.id
FROM house
  LEFT JOIN adresspoint ON adresspoint.id = house.adresspointid
        AND st_contains(st_makeenvelope('562827','6326412','563365','6326590',25832),adresspoint.position)
WHERE adresspoint.id is not null -- ***change is here***
OFFSET 0 LIMIT 100

It runs very slowly (like 100 times slower).

I wonder what is the reason behind it and is there any way I can optimize it? I have indexes on position, ids etc.


texts: Error While Fetching Data From Postgres In Qgis Plugin

texts:
i am trying to fetch data from postgres in qgis and the connection code is working fine in python console but it gives error when i apply that code to my plugin like..

AttributeError: 'NoneType' object has no attribute 'connect'

and my code for connection is

import psycopg2

connection = psycopg2.connect(
    database="test",
    user="postgres",
    password="123",
    host="localhost",
    port=5432
    )

try:
    cursor = connection.cursor()
    check = '''SELECT * FROM testid WHERE "Property_ID" = '54' '''
    cursor.execute(check)
    a = cursor.fetchall()
    print a

except psycopg2.OperationalError:
    pass


texts: How to determine a point lays within a polyline in arc engine?

texts:
I recently asked a similiar question on this, however I needed to know if there is a way to solve this if I only have a IPolyline Feature, and not the whole layer.

I want to know if I can determine if a point lays within my Polyline feature.

I tried using both intersect and difference from ITopologicalOperator, but still no use.

Please Help:

  private bool PointLaysWithinLine(IPolyline currentPath, IPoint iPoint)
    {
        IPolyline newLine = new PolylineClass();
        ((IPointCollection)newLine).AddPoint(iPoint);
        newLine.Project(currentPath.SpatialReference);

//This doesn't work
        IPolyline intersection = ((ITopologicalOperator)currentPath).Intersect(newLine, esriGeometryDimension.esriGeometryNoDimension) as IPolyline;

//This also doesn't work
        IPolyline difference = ((ITopologicalOperator)currentPath).Difference(newLine) as IPolyline;

        return currentPath.Length != difference.Length;
    }

*************************UPDATE*************************

So after playing around for 2 days with user's comments and ideas. I came up with this code (thanks jakub)

point.SpatialReference = IMap.SpatialReference;
ISpatialFilter spatialFilter = new SpatialFilter();
spatialFilter.Geometry = point;
spatialFilter.GeometryField = ((IFeatureLayer)layer).FeatureClass.ShapeFieldName;
spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
IFeatureCursor flowPathCursor = ((IFeatureLayer)layer).FeatureClass.Search(spatialFilter, true);

if (flowPathCursor.NextFeature() != null)
   result = true;

The above code works for MOST time, however, some rare case it still CANNOT find the feature. Note that I USED snapping (vertices/edges) so it should be dead on.

Also I tried using the slope y=mx+b method, and IIdentify2.Identify. They all work MOST time but some rare case (i notice it happening most often on straight EDGES).

Any one have any idea? This is driving me insane.


texts: Algorithm for re-numbering duplicate feature IDs based on distance along line (PyQGIS and QGIS 2.8)

texts:
I've got something working for QGIS 2.8 using Python 2.7 on Windows 7, but I am missing an important part of the puzzle. I'm trying to develop an algorithm to number homes. I've now come to a point where some neighbors have a duplicate house number. For this reason, I've come with a python script that successfully increments the numbers by two (to keep numbers odd or even), depending on the number of duplicates (n-1 for n duplicates, so that exactly one gets to keep the original number). I'm dealing with a metric numbering system, where the house number corresponds with a distance of the home along the road it belongs to.

Starting with this dictionary, where the keys are feature IDs and the values are the concatenated homeNumber_roadNumber:

Deca_dict = {
    "1": "2_506",
    "2": "2_506",
    "3": "2_506",
    "4": "2_600",
    "5": "1_650",
    "6": "2_600"
}

I counted the number of occurrences and kept track of the original number i'm interested in incrementing:

decaAdd_occurrences = {k: (int(k.split('_')[0]), v) for k,v in
                                Counter(Deca_dict.values()).items()}

With the help of a Stack Overflow answer, I adapted my needs and got here:

I adapted the answer's code like so:

changed_deca_dict = {}
alpha_deca_dict = {}
for key, value in Deca_dict.items():
    num, cnt = decaAdd_occurrences[value]
    if cnt > 1:
        route = value.split('_')[1]
        next_num = num + 1
        next_add = str(next_num) + '_' + str(route)
        if not next_add in Deca_dict.values():
            Deca_dict[key] = '{}_{}'.format(next_num, route)
            decaAdd_occurrences[value] = next_num, cnt-1
            changed_deca_dict[key] = '{}_{}_{}'.format(num, next_num, route)
        else:
            alpha_deca_dict[key] = value

This gives me the resulting dictionaries:

Deca_dict -> {
    "1": "3_506",
    "2": "2_506",
    "3": "4_506",
    "4": "3_600",
    "5": "1_650",
    "6": "2_600"
}

changed_deca_dict = {'1':'2_3_506', '3':'2_4_506','4':'2_3_600'}

I have another dictionary with the same keys and their distances along an axis (the order in which I'm trying to number these values)

distance_dict = {'3': 1.0,'1': 2.2,'2': 3.1,'4': 2.0,'5': 1.5,'6': 3.0}

This is the result I would like, re-numbering the duplicate values taking into account the distances of each feature on it's line. Since features 1, 2, and 3 are on the same road, I would re-number the duplicates taking into account their distances along the road. Same for features 4 and 6.

Deca_dict -> {
    "1": "3_506",
    "2": "4_506",
    "3": "2_506",
    "4": "2_600",
    "5": "1_650",
    "6": "3_600"
}

I saw this other thread, which shows how to sort features by a particular field. I'm not sure how to work it in. I'm trying to start by grouping the keys by similar values. I don't want the whole value compared, so this is what I have so far:

for deca_key, deca_value in changed_deca_dict:
    old_num = value.split('_')[0]
    new_num = value.split('_')[1]
    route = value.split('_')[2]
    old_add = str(old_num) + str(route)
    for dist_key, dist_value in distance_dict:
        if dist_key == deca_key:
            check the distance and compare with other Deca_dict.items(), including the one that didn't change
            re-order new_num according to distance along route


texts: Merging overlaying rasters into one using QGIS?

texts:
How do I merge black and white rasters into one file using QGIS?

Also, rasters are overlaying so I´m interested into preserving only black colour.


texts: Overlaying WMS tiles on Google Maps is distorted

texts:
I'm trying to overlay tiles from a WMS server (a map of the ocean biomass) on a Google Map (android), but the projection is slightly distorted.
As you can see in this picture, the longitudes are OK, but the latitudes get distorted the farther away from the equator you get:

When zooming in, the distortion is no longer visible:

The WMS accepts EPSG:4326 bounding boxes (I've reviewed other questions on this topic, but they all use a WMS that accepts different coordinate systems).
I'm calculating the WMS BBOX coordinates out of the (x,y,zoom) parameters in getTile by calling tile2boundingBox(x,y,zoom) (calculation taken from https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames):
double tile2lon(int x, int z) {
  return x / pow(2.0, z) * 360.0 - 180;
}

double tile2lat(int y, int z) {
  double n = pi - (2.0 * pi * y) / pow(2.0, z);
  return vector_math.degrees(atan(sinh(n)));
}

tile2boundingBox(final int x, final int y, final int zoom) async {
  final north = tile2lat(y, zoom);
  final south = tile2lat(y + 1, zoom);
  final west = tile2lon(x, zoom);
  final east = tile2lon(x + 1, zoom);

  return [south, west, north, east];
}

And then I create a WMS url to fetch the tile with the BBOX e.g:
https://my.cmems-du.eu/thredds/wms/global-reanalysis-bio-001-029-daily?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=chl&STYLES=boxfill%2Frainbow&TILED=true&COLORSCALERANGE=0.049000002%2C10.9050001&WIDTH=512&HEIGHT=512&CRS=EPSG%3A4326&BBOX=0.0,-90.0,40.97989806962013,-45.0
I've also compared my calculation to the python program in MapTiler: https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/ and the coordinated are exactly the same:

What am I doing wrong?
Thanks.