texts: GDAL deletes metedata file when creating GTiff

texts:
Reposting question explaining it better.

I am doing TOA corrections on Landsat8 Images in Python(3.5) by reading all image_bands and after processing them creating new GTiff of those images with GDAL, example of code:

#THIS part saves new image
def array2raster(fname,  nx, ny, matriz, geot, proj):
    drv = gdal.GetDriverByName("GTiff")
    dst_ds = drv.Create(fname, nx, ny, 1, gdal.GDT_Byte)
    dst_ds.SetGeoTransform(geot)
    dst_ds.SetProjection(proj)
    dst_ds.GetRasterBand(1).WriteArray(matriz)
    # dst_ds.GetRasterBand(1).SetMetadata({"Band": "%d" % band, "Units": "W/(m2*ster*um)", "Data": "TOA Radiance"})
    dst_ds.GetRasterBand(1).FlushCache()

#Here is made the convertion
def convert_to_radiance(fname, band, gain, bias):
    global img_rad, img_original
    g = gdal.Open(fname)

    gl = g.ReadAsArray()
    print("Dimension: ", gl.shape, gl.dtype)
    radiance = gl * gain[band] + bias[band]

    img_original[:, :, band] = gl
    img_rad[:, :, band] = radiance

    a=fname
    output_fname = a.replace(".TIF", "_TOARAD.TIF")
    print ("Original fname  " , fname)
    print ("Creating output " , output_fname)
    # Create output dataset if `first_time` is true
    nx = g.RasterXSize
    ny = g.RasterYSize
    print("Dimension imagen: ", nx, ny)
    proj = g.GetProjectionRef()
    geoT = g.GetGeoTransform()
    array2raster(output_fname, nx, ny, radiance, geoT, proj)

#this is the main function:

def main():

    global prefix
    metadata_file = get_metadata('D:\Minam\LC80030692016324LGN00\LC80030692016324LGN00_MTL.txt')
    prefix = metadata_file.split("_")[0]
    print ("path de metadata:")
    print (metadata_file)
    (bias_rad, gain_rad, bias_ref, gain_ref,lradmax, lrefmax, nx, ny, sun, desun) = process_metadata(metadata_file) #here i Open the metadata to get values

    start_time = time.time()

    for the_band in [1, 2, 3, 4, 5, 6, 7,  9]:
        input_file = prefix + "_B%d.TIF" % the_band
        print("--------------------")
        print("Procesando imagen: {}".format(the_band))
        print(input_file)

        convert_to_radiance(input_file, the_band, gain_rad, bias_rad)
        convert_to_reflectance(input_file, the_band, gain_ref, bias_ref, sun)

If is the first time processing the bands there is no problem saving the images and the metadata still exists, problem appears when trying to do it again with the images already created, here is when metedata file dissappear, debugging code, it happen on the array2rasterfunction but can't explain why because the metadata file is never called in the function again, what could be causing it?


texts:
to start off with:

Coding in C++

Using WFS Version 1.0.0

Using QT 5.4.1

Windows 8 Latest versions

I'm trying to read in a WFS in code via Gdal compiled with Xerces - I'm able to read the GML files so I can pull spatial point/lines/polygons etc. But I'm not too sure how I do this with the feature names/layers. Any help?]

This is what I'm trying to get: The Cities layer name and Continents layer name, but via code.

<schema xmlns:ms="http://mapserver.gis.umn.edu/mapserver"
xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsd="http://www.w3.org/2001/XMLSchema"
xmlns="http://www.w3.org/2001/XMLSchema" xmlns:gml="http://www.opengis.net/gml"
targetNamespace="http://mapserver.gis.umn.edu/mapserver"elementFormDefault="qualified" version="0.1">
<import namespace="http://www.opengis.net/gml"
schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
<element name="continents" type="ms:continentsType"
substitutionGroup="gml:_Feature"/>

<complexType name="continentsType">
<complexContent>
<extension base="gml:AbstractFeatureType">
<sequence>
<element name="msGeometry" type="gml:GeometryPropertyType" minOccurs="0"maxOccurs="1"/>
<element name="NA2DESC" type="string"/>
<element name="NA3DESC" type="string"/>
</sequence>
</extension>
</complexContent>
</complexType>
<element name="cities" type="ms:citiesType" substitutionGroup="gml:_Feature"/>
<complexType name="citiesType">
<complexContent>
<extension base="gml:AbstractFeatureType">
<sequence>
<element name="msGeometry" type="gml:GeometryPropertyType" minOccurs="0" maxOccurs="1"/>
<element name="POPULATION" type="string"/>
<element name="NAME" type="string"/>
</sequence>
</extension>
</complexContent>
</complexType>
</schema>

You can call the Features type by using this WFS GET request
http://demo.mapserver.org/cgi-bin/wfs?SERVICE=WFS&VERSION=1.0.0&REQUEST=DescribeFeatureType


texts: Check or change default value for domain of attribute in feature class

texts: I want to check (and potentially change) the default value for an attribute domain in a feature class.  However, when I access the feature class properties pop up box and click the Subtypes tab, everything is "greyed out" and I am unable to scroll down to the Field Name I need in the Default Values and Domains box to check the default value.  I am working with an enterprise geodatabase.  Is this a permissions issue?  Or should I be doing something differently?


texts: Getting table of daily temperatures for host of cities using Google Earth Engine?

texts:
I'm not experienced in using Google Earth Engine (GEE), so I've been beating my head against the wall for the past week trying to figure out how to do this. I am trying to export a table with daily average temperatures for cities in 3 select states, from Jan 1, 2018 to Dec 31, 2018. I have been able to get annual averages for each city - but what I really want is a daily average for each city. How should I think about this problem? Should I write a loop? The link to my code is: click here
var gridmet: ImageCollection "GRIDMET: University of Idaho Gridded Surface Meteorological Dataset"

var prism18 = prism.filterDate('2018-01-01','2018-12-31').select('tmean');

var gridmet18 = gridmet.filterDate('2018-01-01','2018-12-31').select('tmmx','tmmn','pr');

var cityMean = gridmet18.reduceRegions(
{
  'reducer': ee.Reducer.mean(),
  'scale': 4000,
  'collection': cities,
});

The asset 'cities' was created in ArcGIS. It contains city boundaries for the states of California, New York and Texas. A lot of the other things that I have tried so far have been commented out.


texts: 1'' DHM: Hillshade is too rough - how to "generalize"?

texts: I got this hillshade from my DHM (I used the hillshade-module from QGIS). For my purpose there is too much relief detail (roughness). How can I get rid of this in QGIS or GRASS?


texts: How to transfer QGIS configurations to latest version?

texts: I have QGIS 2.18.11 on my machine. There are saved wms, wfs and DB connections. After installing QGIS 3.0 I can see saved wms connections but no DB connection is listed in 3.0. Collegue's problem is oppposite, he can list DB connections but cannot list WMS connections. How can we transfer 2.18 settings to 3.0?


texts: Import Point Data to Oracle through FME with considering few other Point file attribute

texts:
I have three Point file (*.mif); Pole, Switch and Transformer

Pole is the main feature which will be my Output, but I need additional information within pole layer attribute using information from Switch and Transformer layer.

Pole layer have 3 attribute field: PoleID, Switch and Transformer
The attribute I will use from Switch layer is SwitchType – (i.e. LV/HV etc.)
From transformer, I want only the status of availability within pole (Yes / No )

Condition:

a. I want to search any Switch which is available within 2mtrs buffer to a Pole

b. If switch identified, then I will read the attribute field of the switch; SwitchType – (i.e. LV/HV etc.) and update in the particular Pole attribute in the Switch field.

c. For Transformer, I want to search any Transformer which is available within 2mtrs buffer to a Pole

d. If Transformer Identified, then the Pole attribute field (Transformer) will updated with Yes otherwise No

e. Finally all the poles with updated attributes will import to Oracle

Kindly help me to update these attribute updation using different point layer


texts: How do I reduce the size of a giant SVG file in QGIS

texts:
I'm using windows 10 X64 and QGIS 2.18.1.  My goal is to create detailed large wall prints using high density contour layers. When I export from print composer as an SVG, I create files from 10 MB to 110MB. I'm using Inkscape to manipulate the stroke/pen width and then subsequently export as a PNG. The resulting lag is disruptive. Since my goals for print size start at ~ 40" X 60" I am unsure of how much detail loss I can handle.

here is an average file (41 MB):
https://drive.google.com/file/d/0B49fTXAI05SweVJqOXJUZ3JTNm8/view?usp=sharing


texts:
I'm trying to calculate Shannon diversity index for trees in a urban context.
I have the following attribute table which I got joining attributes by location between a vector points and a vector grid. The field "res" has the tree species from the point layer. The field "id_malla" has the grid identifier ("id malla" means grid ID in Spanish).

I got the field "act_tot_id" by counting total number of points in every grid cell with the expression:

count("res", "id_malla")

I got the field "act_div_id", which is the number of species by each cell code, by the field calculator with this expression:

count_distinct(“res”, “id_malla”)

The issue comes when I try to calculate the number of individuals within each species by cell code ("id_malla"). I've tried several times with many combinations but I can't get the result. Here are a few of the expressions I've tried:

count("res", "act_tot_id", "id_malla")
count("res", "id_malla", "act_div_id")
count("res", "res", "act_div_id")

None of these gives any proper result

Can someone help with this?


texts: No manifest.json or widget.html file for WAB app

texts: I want wanted to configure a select widget in my WAB code on my local to have only the selectBy Polygon instead of the out of the box SelectBy square settings. So I created a new WAB app with those settings, downloaded the code, and attempted to copy the new select widget from the downloaded folder. I immediately realized there was no Widget.html file and no manifest.json file. Everything else looked similar in terms of file structure. I also noticed that the Widget.js code looked a lot different (It was no formatted correctly and looked almost minified) than the original select Widget.js file. For example, there was a selectBy parameter in the new code, but nothing in the original code (I thought maybe I could just modify the original). Any suggestions as to what's going on here? The new widget definitely won't add to my app view without those missing files.


texts: How to remove a map in QGIS GRASS Plugin

texts:
In QGIS Desktop 2.18.14 with GRASS 7.2.2, I started cleaning vector layers and it worked nicely.

Issue: every time I load a shapefile onto GRASS Plugin by v.in.ogr (v.in.ogr.qgis), it creates a map (its own vector set) within GRASS Database.

How can I remove these maps? Most of them are just temporally files and I do not want them to stay in the harddisk forever.

In earlier documentation there appears a module(?) "GRASS Browser" which could be used to delete maps, or probably g.remove in GRASS (standalone) is the function I need, but I cannot locate them.

EDIT I might have begun with a wrong assumption that we needed to clean mapsets / maps levels. GRASS users can simply remove Location by deleting the folder and start over. It seems fairly practical way to deal with pile of maps.


texts:
I'm working on a qgis project, where I have shapefiles : one is a polygon of the area of interest, the other one is a point type file.

I need to split the polygon file based on a number of points.  For example, every 2000 points will cover one tile of the polygon.
So for 10000 points, I'll divide my polygon into 5 tiles each covering an area of 2000 points.

How can I perform this task automatically, instead of drawing lines to cut polygon everytime ?


texts: Correct long / lat referencing typo in vector layer

texts:
This layer is completely misallinged in terms of the standard WGS84 referencing system. I played around a lot with affine transformations and the like. But I presume it is just a typo in the underlying long/lat (e.g. instead of -655 it should be -6.55).
Is there any way to "retype" long/lat manually in QGIS or gvSIG?


texts: What to filter out if making world map from OpenStreetMap data?

texts:
I'm making a world pin map, high resolution, high detail map.

I've downloaded the planet-latest and played around with it. Before importing the file into a database I need to filter out nodes/objects so that it fits on my harddrive.

I am a bit lost as to what tags can I remove to still be able to make a good world map.
I've been playing around with osmfilter and osmconvert.

Here's an example of how many tags are in the dataset.

  340712803 building
  195127703 source
  145411170 highway
   87781369 addr:housenumber
   81280839 addr:street
   68157172 name
   64141428 addr:city
   58872520 addr:postcode
   35372957 natural
   26512596 addr:country
   25334436 surface
   25106849 source:date
   24669992 landuse
   18837161 power
   16953011 waterway

I'm pretty sure I can remove building, but I have a hard time figuring out what is safe to remove.


texts:
I want to convert cityGML files into shapefiles which show the 3D-objects from above in 2D. It would be also nice to keep the height of the objects in the attribute table or as z-value. I want to use either QGIS, SAGA or GDAL to solve the problem.
If I simply export the GML in QGIS to shapefiles, several polygons are corrupted (when I color them, only the outlines are showing, they are not filled). And I can only keep the height if it was in the GML file added as attribute. If I use ogr2ogr from GDAL, the same problem occurs.
I can repair the polygons with the "fix geometries" tool in QGIS, but I am not sure if this always works or if I would have to check all the polygons everytime to check nothing went wrong...
Here is my current code:
# set path to the gml
gml_path = r"H:\QGIS\original.gml"
intermediate_path = r"H:\QGIS\exported.shp"
shp_path = r"H:\QGIS\repaired.shp"

# load gml
gml_input =  QgsVectorLayer(gml_path, "gml_input", "ogr")

# export as shp
QgsVectorFileWriter.writeAsVectorFormat(gml_input, fileName=intermediate_path, fileEncoding="utf-8", driverName="ESRI Shapefile")

# repair shp
processing.run("native:fixgeometries", {'INPUT':intermediate_path, 'OUTPUT':shp_path})

Are there better or more reliable ways to do this?
Edit: I also tried "cheating" my way around it by trying to clip the GML-objects with a polygon with a bigger extent. This, however, did not work.


texts: PostGIS Spatial query taking a very long time

texts:
I'm trying to return a count of crimes (points) that occur in each police grid (polygons), as well as the geometry data, but the query takes ~22 seconds. This seems oddly high and will be problematic for using it to draw vectors in OpenLayers.

I have 28000 crime points, and 56 grid polygons in my dataset.

This is the query I'm using:

SELECT ST_AsText(ST_Collect(public.cspd_patrol_grid.geom)) AS geo, count(*) AS total_crimes, public.cspd_patrol_grid.command
    FROM public.crimes, public.cspd_patrol_grid
    WHERE ST_Contains(public.cspd_patrol_grid.geom, public.crimes.geo)
    GROUP BY public.cspd_patrol_grid.command;

I thought that it was possibly the size of the geometries that was slowing it down, so I removed that from the SELECT statement, and the query took a more appropriate 2 seconds to return only the total number of crimes in each polygon:

SELECT count(*) AS total_crimes, public.cspd_patrol_grid.command
FROM public.crimes, public.cspd_patrol_grid
WHERE ST_Contains(public.cspd_patrol_grid.geom, public.crimes.geo)
GROUP BY public.cspd_patrol_grid.command;

However, where I'm getting confused is when I do return the geometry, but limit the crimes being queried from all of them to just 'Assault', which makes up about 3% of my crime points. This query only takes 600 ms, even though it still returns 100% of my grid geometries. This makes me think that it has nothing to do with the geometry data, and that there's something wrong with my query.

SELECT ST_AsText(ST_Collect(public.cspd_patrol_grid.geom)) AS geo, count(*) AS total_crimes, public.cspd_patrol_grid.command
FROM public.crimes, public.cspd_patrol_grid
WHERE ST_Contains(public.cspd_patrol_grid.geom, public.crimes.geo) AND public.crimes.crimetype = 'Assault'
GROUP BY public.cspd_patrol_grid.command;

Is there anything noticeably wrong with the original query, or is there something else going on that results in the latter 2 queries being more efficient?


texts:
I am making a map containing information about public transportation: busses, trams, etc. The map will have a layer containing for example the tram lines. Each line has its own colour, and is represented by a line string geometry.

The problem is, many of the lines have overlapping parts, where more than one tram line cover the same section. To show this to the user, I'd rather want the lines to run parallel to each other instead of being drawn on top of each other. For an example, see how Google Maps show the New York subway lines.

I suspect this is a quite common problem in cartography, but don't know what terminology I should search for.

I am using PostGIS/GeoServer/OpenLayers as my stack, but any open source solution would be acceptable.


texts: Dissolving/merging polygons taking into account area of intersection

texts:
How in QGIS do I dissolve/merge polygons taking into account the area of intersection?

For example, if area of intersection less then 15% of polygon area then operation should't apply.


texts: Combine data and geometry from vector layers without merging

texts:
I have 3 line type vector layers representing municipal infrastructure; sanitary, storm, and water main pipes.  Each of these layers should remain separate and all of these layers have different fields, but some common attributes.
I'm trying to create a single atlas that looks at a few of the common attributes between the 3 layers to highlight missing information to present to the client; I don't want to merge the tables to do so or create 3 separate atlases.
I've created a virtual layer that selects the proper information from each table and provides the results in a single table - but unless there's something I'm missing, I cannot include any geometry in the virtual layer, and therefore my atlas doesn't work.


texts: Unable to Configure Pop-Ups in ArcGIS Online

texts:
A particular map of mine published to arcgis portal is not allowing me to configure pop-ups. I can enable the pop-ups but when i click 'configure' nothing happens what-so-ever.

Anyone have any idea what could be causing this?


texts:
How to correctly select a edge of a polygon in Google Maps?

So far, I was able to create a simple demo here

You can try the code to see that there's a bug wherein it doesn't select the last edge of the polygon.
Any idea how to deal with this correctly?

var markerDist;
var startPointMarker = new google.maps.Marker();
var endPointMarker = new google.maps.Marker();
var map;

var clickedPolylineLatLng = new google.maps.MVCArray();
var clickedPolylineArray = [];
var clickedPolyline = new google.maps.Polyline({
  path: clickedPolylineLatLng,
  strokeColor: 'blue',
  strokeWidth: 10
});

function distanceBetween(point1, point2) {
  return google.maps.geometry.spherical.computeDistanceBetween(point1, point2);
}

function initialize() {
  var myLatLng = new google.maps.LatLng(0, -180);
  var myOptions = {
    zoom: 2,
    center: myLatLng,
    mapTypeId: google.maps.MapTypeId.TERRAIN
  };

  map = new google.maps.Map(document.getElementById("map_canvas"), myOptions);

  var flightPlanCoordinates = [new google.maps.LatLng(-27.46758, 153.027892), new google.maps.LatLng(37.772323, -122.214897), new google.maps.LatLng(29.46758, 88.027892), new google.maps.LatLng(20.46758, 97.027892), new google.maps.LatLng(17.772323, 78.214897)];
  var flightPath = new google.maps.Polygon({
    paths: flightPlanCoordinates,
    strokeColor: "#FF0000",
    strokeOpacity: 1.0,
    strokeWeight: 10
  });

  google.maps.event.addListener(flightPath, 'click', function(el) {
    markerDist = {
      p1: '',
      p2: '',
      d: -1
    };

    startPointMarker.setMap(null);
    endPointMarker.setMap(null);

    clickedPolyline.setMap(null);
    //clickedPolylineLatLng.clear();

    var latLng = el.latLng;

    console.log('getPath array');
    console.table(this.getPath().getArray());
    console.log('getPaths array');
    console.log(this.getPath().getArray().length);
    console.log(this.getPath().getLength());

    var allPoints = this.getPath().getArray();
    for (var i = 0; i < allPoints.length - 1; i++) {
      var ab = distanceBetween(allPoints[i], latLng);
      var bc = distanceBetween(latLng, allPoints[i + 1]);
      var ac = distanceBetween(allPoints[i], allPoints[i + 1]);

      if ((parseFloat(markerDist.d) == -1) || parseFloat(markerDist.d) > parseFloat(Math.abs(ab + bc - ac))) {
        markerDist.p1 = allPoints[i];
        markerDist.p2 = allPoints[i + 1];
        markerDist.d = Math.abs(ab + bc - ac);

        clickedPolylineLatLng.clear();
        clickedPolylineLatLng.push(allPoints[i]);
        clickedPolylineLatLng.push(allPoints[i + 1]);

      }
    }
    startPointMarker.setPosition(markerDist.p1);
    endPointMarker.setPosition(markerDist.p2);
    startPointMarker.setMap(map);
    endPointMarker.setMap(map);
    clickedPolyline.setMap(map);
  });
  flightPath.setMap(map);
}

window.onload = initialize;


texts: Populating polygon extents to attribute table in QGIS?

texts:
Is there a way to calculate the XMIN, YMIN, XMAX and YMAX of a polygon and write it to a field in an attribute table?

The Geometry functions available in the Field Calculator are for points.

I am using 1.8.0 for Windows.


texts: How to have newly added lines split existing lines where they intersect, using QGIS

texts:
I have a line shapefile where the lines form a closed polygon.  Using QGIS (2.8 on OSX 10.7) I would like to add lines to the layer, basically splitting the polygon into quadrants.  Then, I will need to convert the line shapefile to a polygon shapefile, so that the newly formed quadrants will be separate polygons.

I am encountering some difficulties, and am hoping I can be pointed in the right direction.  First, when I add the new lines, they don't "split" the existing boundary lines of the polygon where the new lines intersect (I am snapping the new lines to existing vertices of the boundary lines).  Second, where the new lines cross in the middle, they also don't split/intersect.

In order to create a polygon shapefile with the necessary quadrants, it seems that the lines need to be intersected, or else I only get the main original polygon as a result (using the Polygonize process).  I have tried the "Split lines with lines" option (using my line shapefile as both the input and output layer), and I have tried v.clean from the GRASS processes, but neither of these approaches split the lines where they intersect (either in the middle, or on the boundary) in all cases, only some.

I feel that the easiest way would be if the new lines would split the existing lines when I add them (and also split each other where they cross), but I can't figure out how to do that.


texts: Examples of thematic mapping web services built with FOSS GIS

texts: What are some current examples of web services serving thematic mapping based on free and open source software?  The thematic mapping could be in the form of choropleth or point layers.  It would be very interesting to see some good publicly accessible geo-information systems, with the emphasis as much on simple but well-implemented projects as on those that are technically accomplished.


texts: GetFeatureInfo with leaflet.extras2 and Shiny

texts:
I am building an R Shiny App that utilizes the following WMS. However, the GetFeatureInfo capabilities for click-enabled popup (leaflet.extras2 functionality) returns a summary of the GetCapabilities. Changing the GetCapabilities to GetFeatureInfo request for the WMS url does not provide any functionality. I have contacted the WMS admin and they GetFeatureInfo is reportedly enabled. Any reasons why it wouldn't be working in with this WMS? I have no previous experience working with WMS.
Additionally, the WMS layer does not show in the map window when this popup (image attached) functionality is enabled. However, ending the WMS url at: https://ortho.gis.iastate.edu/arcgis/services/ortho/lidar_dem_3m/ImageServer/WMSServer allows map visualization but not popup functionality.
library(leaflet)
library(leaflet.extras)
library(leaflet.extras2) # built in GetFeatureInfo
library(htmltools)
library(htmlwidgets)
library(XML)

leaflet() %>%
  setView(lng = -93.6319,
          lat = 42.0308,
          zoom = 12) %>%
  leaflet.extras2::addWMS(
    # "https://ortho.gis.iastate.edu/arcgis/services/ortho/lidar_dem_3m/ImageServer/WMSServer", # working alternative
    "https://ortho.gis.iastate.edu/arcgis/services/ortho/lidar_dem_3m/ImageServer/WMSServer?request=GetCapabilities&service",
    layers = "lidar_dem_3m",
    options = WMSTileOptions(
      format = "image/png",
      version = "1.3.0",
      transparent = T,
      info_format = "text/xml"
    ),
    group = "dem3",
    popupOptions = popupOptions(maxWidth = 300, closeOnClick = T)
  ) %>%
  addLayersControl(baseGroups = c("dem3"))


texts: ClassifyPixelsUsingDeepLearning Output Error from ArcGIS Pro

texts:
While producing the final pixel block, which is a dictionary of arrays of predictions, ArcGIS Pro errors with

Unable to unpack the python NumPy array associated with the specified
raster into a pixel block. [output] The shape of the specified array
is unexpected. [pixels: incoming 1 x 32 x 32, expected: 4 x 32 x 32.]

However, the true expected output is a 1x32x32, as it produces a class for each pixel in the image tile.


texts: Updating the reference data for an existing geocode service

texts:
I created an address locator in an ArcSDE gdb, based on data stored in a point feature class in the same SDE gdb. I then shared the locator as a geocode service and am consuming it in a Silverlight app. All of this works perfectly.

I have now updated my reference data in the point feature class, and rebuilt the locator. Would I have to recreate the geocode service to access the updated reference data? I have been reading the extensive help on registering data with ArcGIS Server, and followed the tutorial on publishing a geocode service, but it appears that the locator data has been copied over to the server, even though the SDE gdb is accessible to it.

I tried replacing the service, but it wouldn't let me (had to specify a new name). When I then tried connect to the second geocoding service, it said connection failed (even though it is set up the same way as the original service).


texts: ESRI Feature Service - Unable to complete operation

texts:
I have a road network published on ArcGIS Portal 10.6. I created it using ArcMap 10.6.

While publishing I enabled the network service and later I changed the service properties to include the feature access and allowed the update operation.

When I try updating attributes of a feature in roads layer it is not updating. I checked the network tab in Developers tool and i noticed the update operation is giving the following response

{"error":{"code":400,"message":"Unable to complete operation.","details":["Internal error during object update."]}}

Can I edit a layer once it is part of road network and is published?

Update:

I checked the ArcGIS Server logs after trying the edit function and it gave me the following error. I only edited roads layer but looks like it is trying to edit barriers as well.

SEVERE  An error occurred.  RoadNetwork.MapServer
SEVERE  Error: Insufficient permissions [ND_14394_DirtyAreas][STATE_ID = 393].  RoadNetwork.MapServer
SEVERE  Initialization of Layer: Polygon Barriers failed.   RoadNetwork.MapServer
SEVERE  Layer: Polygon Barriers does not use SDE GeoDatabase. Feature Access capabilty requires the data to be on SDE GeoDatabase.  RoadNetwork.MapServer
SEVERE  Initialization of Layer: Line Barriers failed.  RoadNetwork.MapServer
SEVERE  Layer: Line Barriers does not use SDE GeoDatabase. Feature Access capabilty requires the data to be on SDE GeoDatabase. RoadNetwork.MapServer
SEVERE  Initialization of Layer: Routes failed. RoadNetwork.MapServer
SEVERE  Layer: Routes does not use SDE GeoDatabase. Feature Access capabilty requires the data to be on SDE GeoDatabase.    RoadNetwork.MapServer
SEVERE  Initialization of Layer: Point Barriers failed. RoadNetwork.MapServer
SEVERE  Layer: Point Barriers does not use SDE GeoDatabase. Feature Access capabilty requires the data to be on SDE GeoDatabase.    RoadNetwork.MapServer
SEVERE  Initialization of Layer: Stops failed.  RoadNetwork.MapServer
SEVERE  Layer: Stops does not use SDE GeoDatabase. Feature Access capabilty requires the data to be on SDE GeoDatabase. RoadNetwork.MapServer
SEVERE  This exception was thrown after the response was committed. Access to this resource is not allowed  Rest
SEVERE  Response already committed. Cannot forward to error page.   Rest
SEVERE  Response already committed. Cannot forward to error page.   Rest


texts: Field Calculator Parsing error SyntaxError: invalid syntax (line 2)

texts:
I'm writing a conditional in the Field Calculator and it's giving me an error of syntax in line 2.
Can someone see it?

The Code Block is this:

def ZFactor(LAT_CENTR) :
  if (10.000 < math.fabs(LAT_CENTR) and math.fabs(LAT_CENTR) => 00.000) :
    return 0.00000898
  elif (math.fabs(LAT_CENTR) => 10.000 and math.fabs(LAT_CENTR) < 20.000):
    return 0.00000912
  elif (math.fabs(LAT_CENTR) => 20.000 and math.fabs(LAT_CENTR) < 30.000):
   return 0.00000956
  elif (math.fabs(LAT_CENTR) => 30.000 and math.fabs(LAT_CENTR) < 40.000):
     return 0.00001036
  elif (math.fabs(LAT_CENTR) => 40.000 and math.fabs(LAT_CENTR) < 50.000):
    return 0.00001171
  elif (math.fabs(LAT_CENTR) => 50.000 and math.fabs(LAT_CENTR) < 60.000):
    return 0.00001395
  elif (math.fabs(LAT_CENTR) => 60.000 and math.fabs(LAT_CENTR) < 70.000):
    return 0.00001792
  elif (math.fabs(LAT_CENTR) => 70.000 and math.fabs(LAT_CENTR) < 80.000):
    return 0.00002619
  elif (math.fabs(LAT_CENTR) => 80.000 and math.fabs(LAT_CENTR) <= 90.000):
    return 0.00005156
  else:
    return 1


texts: PostGIS point distance from polygon not 0 when on edge

texts:
I have a pretty simple problem:

SELECT * FROM ST_Distance(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -39)'), False)

returns 118.45656881
Shouldn't it be 0? Isn't the point on the edge of the polygon?
Why does

SELECT * FROM ST_Distance(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -40)'), False)

return 0?

SELECT * FROM ST_Distance(ST_GeographyFromText('POINT(70 -39.5)'), ST_GeographyFromText('POINT(70 -39)'));
SELECT * FROM ST_Distance(ST_GeographyFromText('POINT(70 -39.5)'), ST_GeographyFromText('POINT(70 -40)'));

Why does this return two different (although only slightly) numbers?


texts:
Trying to import a medium-sized .csv file (ten columns, approximately 500 lines) into QGIS.

I have a database in Microsoft Access, and am exporting files from that as .csv to import into QGIS. I have seen the problems using Excel, and so am doing it via OpenOffice. I have also tried making them into .tsv files via Google Sheets.

When I press 'add delimited text layer', the preview window truncates my file to 20 lines. When I then select the fields for latitude and longitude (I'm using UTM), it imports it truncated to (a different) 15 lines. It also changes the latitude and longitude values to either zero or NULL.

The file looks fine in Notepad. If I make it import it without geometry, it still truncates it in the preview window, but then imports the whole table. However, I then can't put the points on the map!

Does anyone know what is going on, or any way to work around this other than creating a new shapefile and manually placing every single point?

EDIT: The data looks like this:

ID,Region,Site,Type,Number,Nickname,Secure Dates,Uncertain Dates,Latitude,Longitude,Notes
1,Messenia,Koryphasio,Tholos,0,,"LH I, LH IIA",MH,"4,094,405.15","560,361.67",
5,Messenia,Pylos,Tholos,5,Grave Circle/Vagenas,"LH I, LH IIA, LH IIB, LH IIIA1",MH,"4,097,946.59","561,721.62",
6,Messenia,Voidokoilia,Tholos,0,,"LH I, LH IIA, LH IIB, LH IIIB1, LH IIIB2",,"4,091,240.82","558,773.89",

And so on all the way down.


texts: Export basemap and vector layers into a single image with OpenLayers 2.11?

texts:
Has anyone been able to export OpenLayers 2.11 map and features as a static image?

This example is exactly what I'd like to do, but it's using OL3 alpha which is not an option for me...

http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/examples/exportMapCanvas.html

Thanks!
Vincent


texts: can QGIS for Mac OSX 10.6 handle (both read and generate) GeoPDF files?

texts: I would like to know if QGIS for Mac OSX 10.6 can handle (both read and generate) GeoPDF files. Can anybody help ?


texts: Cut geoprocessing not outputing expected results

texts:
I wish to cut a vector layer with another layer (single polygon).
I have tried both QGIS and GDAL cut geoprocessing with little satisfaction.

Input layer:

Output layer (with cutting polyon in thick black border):

See the blank space that should be green.

Same operation using ArcGIS gives correct result. (Sorry I cannot paste screenshot here cause I lack reputation... you'll have to believe me)

Output message says one or more geometry invalid
or TopologyException: Input geom 0 is invalid: Self-intersection at or near point.

Any idea how to tackle this?


texts: Using YSLD functions and creating legends

texts:
My goal is to create legends, for WMS layers, which can be shown in a JavaScript/OpenLayers 3 client.

Furthermore I found out, that GeoServer can generate legends via the SLD styles. This legends can be requested with a GetLegendGraphic-call.

I chose YSLD over SLD to write my styles because it is more readable and easier to write. Therefore my YSLD files look more or less like this:

name: simple-poly
title: 'A simple ysld style'
feature-styles:
- rules:
  - name: Rule uses Categorize()
    scale: [min,max]
    symbolizers:
    - polygon:
        fill-color: ${Categorize(foo,
          '#006600','73',
          '#00CC00','79',
          '#66FF66','85',
          '#996600','88',
          '#FFCC66','94',
          '#CCCCCC')}

Especially the Categorize() function saves a lot of lines. With these few lines my WMS layer gets rendered as expected.

The drawback with this approach seems to be that GeoServer can not handle the Categorize() function when creating a legend. Because the result of the GetLegendGraphic-call is an image with a gray square:

What is the best way to create a legend using YSLD, especially when Categorize() is used? Or is it simply not possible to use Categorize() and creating legends automatically?


texts:
I am trying to get MongoDB working with GeoServer. I am using the BoundlessGeo Suite (comes with GeoServer and other tools) I was able to find the MongoDB extension for GeoServer. Got it installed according to these instructions:

Getting MongoDB Working with GeoServer

I found this version of the extension:
http://repo.opengeo.org/org/geotools/gt-mongodb/12-SNAPSHOT/

I am able to see MongoDB as an option in the Data Store Menu. MongoDB is also installed on the server locally for testing. I was able to create a database without much issue.

After putting in the Data Store Settings, I attempt to upload a shapefile using GeoExplorer when I get the following error:

"Name is required for PropertyType"

When I searched for the error I found almost no results but what I did find seemed to indicate that I needed to make an XML file to designate where the Mongo Info is. I would just have no idea about what needs to go in that file or where it needs to go.

http://docs.geoserver.org/stable/en/user/data/app-schema/mapping-file.html

Kind of stuck and wondering if anybody else out there was able to get it to work on any of their setup.


texts: Exporting true color sentinel image in Google Earth Engine?

texts:
I tried to download a sentinel image, but it seems not to work. It downloads a one band image and not 4,3,2 rgb combination as coded.

This is the code I'm using that just downloads 1 band:

var marco = /* color: #d63000 */ee.Geometry.Polygon(
    [[[-70.33198356628418, -16.667604212318242],
      [-70.33138275146484, -16.64129056542051],
      [-70.3795337677002, -16.641043857880046],
      [-70.3795337677002, -16.668426455517558]]]);

var sentinel = ee.Image('COPERNICUS/S2/20180726T144731_20180726T145459_T19KCB').select(['B4', 'B3', 'B2']);

Export.image.toDrive({
 image: sentinel,
 description: 'imageToDriveExample',
 scale: 10,
 region: marco
});


texts: Filter ArcSDE layers of the TOC into a combobox

texts:
I've a list of layers of the current workspace in a combobox, but I want to list only the layers of the ArcSDE Connection, not the others I've in the TOC.
How can I filter the name of this layers?


texts: QGIS (3.10.7) SAGA tools resulting in "The following layers were not correctly generated."

texts:
I am using QGIS 3.10.7 and want to use SAGA tools. I am using an ASCII file (.asc) with a projected coordinate system - EPSG:26919 - NAD83 / UTM zone 19N - Projected.
When attempting to use the tools, they always result in a message like this -

"The following layers were not correctly
generated.C:/Users/Admin/AppData/Local/Temp/processing_7cf0f5b0db874bb9a6b8fe84b44cb047/2a6b2113845f4739b4dec0f12c451504/C_PROF.sdat"

I read a similar post, that recommended "[save to temporary file]." I tried that and it still does not work. The following layers were not correctly generated. Grid. I also don't believe that I have any spaces in my file extension as seen in the error message.
Despite the error, it ends up outputting values. The values seem reasonable, but I am of course suspicious of them given that it says they were not correctly generated.
I first had this issue on my Mac. I tried multiple versions of QGIS and still had the issue. Now I am using a PC with 3.10.7 - still getting the same message.
Any ideas?


texts: QGIS cannot import SciPy (QGIS 3.6 Noosa, Ubuntu 18.04, EnMapBox plugin)

texts:
I had been struggling to get QGIS installed because (currently) the ubuntugis repo seems to not be working well (couldn’t load Python at all). I got around that problem by installing QGIS using flatpak.

I need to use the ENMAPBox plugin, which requires SciPy, Numpy, and Sklearn modules. Ubuntu came with Python 3.6.8 installed and QGIS brought with it Python 3.7.3. ENMAPBox cannot start because it can’t find either of those modules, and importing them from the console yields the same response. I followed the plugin’s instructions to install those packages using python3 -m pip install.

I’ve verified that the executable used by QGIS python console is /usr/bin/python3, and inspecting sys.path shows that /usr/lib/python3/dis-packages is not in the path (so it would make sense that it can’t find the packages I installed using pip3).

sys.path looks like this:

['/app/share/qgis/python',
'/home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python',
'/home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins',
'/app/share/qgis/python/plugins',
'/usr/lib/python37.zip',
'/usr/lib/python3.7',
'/usr/lib/python3.7/lib-dynload',
'/usr/lib/python3.7/site-packages',
'/app/lib/python3.7/site-packages',
'/app/lib/python3.7/site-packages/numpy-1.14.3-py3.7-linux-x86_64.egg',
'/app/lib/python3.7/site-packages/MarkupSafe-1.0-py3.7-linux-x86_64.egg',
'/home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python',
'/home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages',
'/home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapbox',
'/usr/lib/python3/dist-packages']

Note the last item is there because I appended it, it’s not there by default.

To get QGIS to import SciPy I’ve tried a number of things:

Appending /usr/lib/python3/dist-packages to sys.path
Redirecting pip3 to /usr/lib/python3.7/site-packages
Adding .pth files to /usr/lib/python3.7/site-packages and /home/nav00a/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages with links to /usr/lib/python3/dust-packages and /usr/lib/python3.7/site-packages

Unfortunately the behaviour remains the same. If I import SciPy this is what I get:

Traceback (most recent call last):
  File "/usr/lib/python3.7/code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/app/share/qgis/python/qgis/utils.py", line 672, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'scipy'

Am I missing something that’s very obvious?


texts:
I'm using arcgis 10.0. I red the ESRI doc to setting up script tools parameters, but I don't succeed at setting them up.

My parameters before (without variable)

#arcpy.env.workspace = r'Z:\Documents\SIG\tests\synthese_tocorrect.gdb'
#pochoir = r'Z:\Documents\SIG\tests\decoupe.shp' #the clip layer with multiple polygons
#outFolder = r'Z:\Documents\SIG\tests\sorties'

My parameters after, with issues :

import arcpy, os
from arcpy import env
inWorkspace = arcpy.GetParamterAsText(0)#Refers to the geodatabase
pochoir = arcpy.GetParamterAsText(1) #Refers to the existing clip layer
outFolder = arcpy.GetParamterAsText(2) #Output folder.
env.workspace = inWorkspace
arcpy.MakeFeatureLayer_management(pochoir, "pochoirlyr") #Create layer to allow for selection
fcList = arcpy.ListFeatureClasses() #the feature class to be clipped by each of the polygon of the clip layer
OIDField = arcpy.Describe("pochoirlyr").OIDFieldName #Get name of OID field

for fc in fcList:
    arcpy.MakeFeatureLayer_management(os.path.join(arcpy.env.workspace, fc), "lyr")
    pocursor = arcpy.SearchCursor("pochoirlyr")
    for porow in pocursor:
        sql = '"' + OIDField + '" = ' + str(porow.getValue(OIDField)) #SQL to select one feature
        arcpy.SelectLayerByAttribute_management("pochoirlyr", "", sql) #Select polygon feature by OID
        arcpy.SelectLayerByLocation_management("lyr","INTERSECT","pochoirlyr") #Select feature layer to be clipped by selected clip feature
        if arcpy.Describe("lyr").FIDSet: #Check for selection
            out_poly = os.path.join(outFolder, fc + "_" + str(porow.getValue(OIDField))) #Output feature class
            arcpy.Clip_analysis(fc, "pochoirlyr", out_poly)
    del porow
    del pocursor
    arcpy.SelectLayerByAttribute_management("pochoirlyr", "CLEAR_SELECTION")
    arcpy.Delete_management("lyr")

And in the properties of the script tools :

geodatabase : Workspace or Feature Dataset, input, no multivalue...
pochoir : shapefile, input, no multivalue...
sorties : Folder, output...

But when I choose a different parameter when opening the script tool GUI, it does not work. It does the job with the IDE parameters.

Any help ?


texts:
I have docker installation of GeoNode with all defaults: my_geonode, my_geonode_data etc. Now I'm truing to change default password for user my_geonode in PostgreSQL and I did. Also, I have change in .env file to new password
GEONODE_DATABASE=my_geonode
GEONODE_DATABASE_PASSWORD="my new password"

I'm not sure if this is an issue but after changing password geonode is not able to start:
psql: FATAL: password authentication failed for user "my_geonode"
my_geonode is unavailable - sleeping

I'm not sure if this is an issue or not but I do not know where else this new password should also be placed.


texts:
When I attempt to “Export as PDF” a Print Composer image, it hangs up until I press cntrl/alt/del, whereupon I learn that QGIS is not responding and needs to be closed. I’ve tried this as jpg, jpeg, and tif.

Don't know if this is related or a different issue--When I attempt to “Export as PDF”, it goes better until I open the new image in Acrobat Reader. Then I get a message that “An error exists on this page. Acrobat may not display the image correctly.” And indeed, all I get is the background raster (aerial photo) that my map is based on—none of the vectors, etc., that I have overlaid. I've seen in this forum that labels can be an issue, and to use the "new" labels from the ABC icon. Does this refer to labels created in Print Composer (which I have many of), or labels attached to layers (which I don't have any of)? I don't see any ABC icon, and don't understand what this means and cannot find further information.


texts: Mosaic huge numbers of small rasters using ArcGIS

texts:
I have around 20'000+ TIFF-rasters that I need to somehow bring together into one raster. The rasters are model outputs and usually aren't big in size (<< 1 MB each). They do not overlap. The cell size for all of them is 10 m. I have written a little python script that copies the rasters, renames them and then starts the mosaic to new raster process (code snippet below).

######## 4. Umbenennen und kopieren der Outputraster, inkl. löschen der nicht mehr benötigten Verzeichnisse
arcpy.AddMessage("Umbenenn- und Kopiervorgang wird gestartet...")
for raster in listRasters:
    new_name_raster = Sammelordner + raster[-25:-17] + ".tif"
    arcpy.CopyRaster_management(raster, new_name_raster)
    arcpy.Delete_management(raster[:-17])

arcpy.AddMessage("Umbenenn- und Kopiervorgang abgeschlossen.")
kachelliste = arcpy.ListFiles("kachel*")
for item in kachelliste:
    arcpy.Delete_management(item)
arcpy.AddMessage("Leere Verzeichnisse gelöscht.")

######## 5. Mosaikieren der kopierten und umbenennten Raster:
if check_mosaic.lower() == 'true':
    arcpy.AddMessage("Mosaikiervorgang wird gestartet...")
    arcpy.env.workspace = Sammelordner
    arcpy.env.scratchworkspace = Sammelordner
    arcpy.env.overwriteOutput = True

    mosaic_candidates = arcpy.ListRasters("*.tif","TIF")
    mosaic_name = desired_outputname + ".tif"
    arcpy.MosaicToNewRaster_management(mosaic_candidates, Mosaikordner, mosaic_name,"","32_BIT_FLOAT", 10, 1, "MAXIMUM", "LAST")
else:
    arcpy.AddMessage("Mosaikiervorgang nicht durchgeführt.")

The script is now running since more than a week, the mosaic is slowly becoming larger and is about to reach it's goal size (I've tested the script with a smaller amount of rasters and then upscaled the resulting mosaic in order to estimate the goal size of the 20'000+ rasters). So the approach is working but I'd like to speed it up, since there is a chance I'll have to repeat the simulations and then create another mosaic...

The following ideas have been around:

Would creating a mosaic dataset or a raster catalogue be much faster? I thought they weren't the right choices for my problem... I've read this (Fast way to merge 8000 raster?) but didn't find the answers too helpful for my problem.
Would creating smaller mosaics (say i.e. 200 x 100 rasters) and then mosaic these togehter be much faster?
I have read about the use of GDAL, but have no clue about how to access it using ArcGIS (I don't have QGIS or any other opensource GIS applications available at work). Would this help?
Any other suggestions?


texts: .leaflet-control-layers-expanded style control appearance

texts:
This might be kind of a nothingburger, but it's nagging enough that I can't fix it and I'd like to know why this is happening!!!  I'm also wondering if this little bug might be a sign something isn't quite right with my js.

Coding from this example, my layer control doesn't look like it does in the example map.  The only difference I have between the code is using the topographic map instead of the grayscale.  Otherwise my code is identical.  The appearance is the same in both FF and Chrome.

In fact, just to get it to look like this, I had to overwrite the leaflet.css with this:

<style>
    .leaflet-control-layers-expanded {
        width: 125px !important;
    }
</style>

Without this CSS overwrite, the expanded layer control had all 3 items on the same line, under the horizontal rule.

So, how do I get the 2 radio buttons that control the basemap toggle above the horizontal rule in the expanded layer control?


texts: Google maps dynamically load marker after center_changed with minimum distance

texts:
I would like to load dynamically markers (trees) to get displayed around the center of a google map. I call the script xml_trees.php on our server passing the southwest and the nordeast coordinates of the map bounds to get only the trees inside this bounding box. This works fine.

When I'm outside in the field with a tablet displaying my actual gps position and updating the viewport according the gps position, there are to many calls to the database, since even minimal position changes fire the event.
Is there a way to reduce the calls so that the call of the xml_trees.php would only be called, when a changed position exceeds a certain distance (e.g. 10m)?

function getTreesAroundMe(){
    google.maps.event.addListener(map, 'center_changed', function() { // instead of 'bounds_changed'
                  if (map.getZoom()>18)  {
                      var bounds =  map.getBounds();
                      var ne = bounds.getNorthEast();
                      var sw = bounds.getSouthWest();
                      dlurl = "./mapdata/xml_trees.php?minx=" + sw.lng() + "&miny="+sw.lat()+"&maxx=" + ne.lng()+"&maxy="+ne.lat();
              getTrees(dlurl,true);
                  }
    });
}


texts: How to avoid exceeding latitude and longitude limits, when importing a .csv file?

texts:
I am currently working on a project with "on the fly"-projection (WGS 84 / UTM zone 32U).

I want to import a csv. file with points that have a longitude and latitude using the "Add Delimited Text Layer"-tool. In the tool I have to select a CRS. I tried this with several CRS's.

At first, the points will not show up correctly, but if I manually select a random CRS, then the points will show up correctly.

If I zoom out, the following Error will appear:

Error: "latitude or longitude exceeded limits"

I already read, that this might happen because of the lon&lat and a projected CRS, but this does not help me up to now. I think that I might have a lack of knowledge with determining what CRS to assign'.

edit: This is how my .csv file looks like:


texts: Can I set individual opacity levels for an ArcGISDynamicMapServiceLayer?

texts:
I've added a ArcGISDynamicMapServiceLayer and I've used the imageParameters.layerIds property to only show certain layers.

Is there a way to set a different opacity level for each layer?


texts: How to build a function that extracts information from a set of polygons at multiple points in R

texts:
I am new to R, and trying to build a simple function that calculates an average (mean) of Sum_MVKT extracted from a set of predefined point locations.

For this, I have two shapefiles named Traffic_Val (polygon), and StnLocations (points).

"Traffic_Val" contains a large number of polygons and each polygon has an unique value (Sum_MVKT).

"StnLocations" consists of multiple points that lie on top of polygons in "Traffic_Val".

Provided the predefined StnLocations, I'd like to extract the value from Traffic_Val at each station location and calculate the average.

Here is the R script that I have written so far..

library(sp)
library(maptools)

setwd("C:/Users/calculateTraffic")

# load data:
SA <- readShapePoly('Traffic_Val.shp') # This file contains a list of traffic values named as "Sum_MVKT"
proj4string(SA) <- CRS("+init=epsg:26915")

### Existing Station Locations ###
PtsLocation <- readShapePoints('stnLocations.shp')
names(PtsLocation@data)
names(PtsLocation@data)<-c("ID", "x", "y", "Site_desc")

plot(SA)
points(PtsLocation, pch = 1, col = 3)

###########################
#### calculateTraffic #####
###########################

calculateTraffic = function (PtsLocation)
{
  obs = PtsLocation
  if (missing(formulaString) || is.null(formulaString)) {
    eq = dum ~ 1
  } else eq = formulaString
  ############ changed following 'if statement'
  if (!"data.frame" %in% getSlots(class(obs)) & (terms(eq)[[3]] ==
                                               1 || all(all.vars(eq)[-1] %in% dimnames(coordinates(obs))[[2]]))) {
    obs = SpatialPointsDataFrame(obs, data = data.frame(dum = rep(1,
                                                              dim(coordinates(obs))[1])))
    names(obs) = as.character(eq[[2]])
  }

  meanTraffic <<- mean( )

  return(meanTraffic)

}


texts: OpenLayers 3 - in-house tiled ArcGIS basemap not visible

texts:
I am new to openLayers.  I am trying to display an in-house cached basemap from an  ArcGIS REST service.
I can see the basemap if I use "ol.source.TileArcGISRest" for the source property, but the performance is very slow and the map is rendered very blurry.

I am trying instead to use "ol.source.XYZ" for the source property of the layer (as recommended in OpenLayers docs), but I am now not seeing any map.  The X and Y values I see in the Network tab of the developer tools are not in the range given in the REST endpoint for the tiled map.
Our in-house maps all use NY State Plane (EPSG 2263).

FOr instance, at a zoom level of 0 (in the 'view' property' of the code), my sample attempts to get a tile with Z Y X parameters 9/192/150 and 9/192/151, but the start and end Z Y X values at zoom level 9 (as set in this tiled map REST service) are 9/646969/537582 and 9/657961/548957.

I am guessing this is a projection problem, but I thought I had set the default projection in my code. The full example is below:

Both the projection extent and extent in the layer properties are currently set from the extent given in the REST endpoint.  I tried swapping out the projection extent with that given in the epsg URL, and tried commenting either of them out, to no effect.

<!DOCTYPE html>

    <title>Tiled ArcGIS MapServer</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.13.0/css/ol.css" type="text/css">

     script src="http://cdnjs.cloudflare.com/ajax/libs/proj4js/2.3.6/proj4.js" type="text/javascript"></script>


texts: Create DEM with 1 foot contours from lidar data with 1.4 meter resolution

texts:
I have Lidar data which contains multiple .las files for an entire quad. I have figured out which .las is the area I am looking for. I was told the data is in 1.4 meter resolution and would work to create a contour map with 1 foot intervals. So that's what I have to work with.

I am not sure what the 1.4 meter resolution means or how it helps to create 1 foot intervals. Would the first step be to use the las2dem tool? And what would I set the step size/pixel size to?

I have created a DEM by inputting my .las file, no filter, pixel size:1. I then used raster>extraction>contours and input my DEM with interval set to 1. Is this correct for 1' intervals, or is there something I need to do to compensate for the 1.4 meter resolution of the .las?

I am mostly confused on how to make sure the intervals are set to 1 foot and not 1 meter, and what the 1.4 meter resolution means.

Also, I am working in QGIS and using the LAStools set.


texts: How to remove overlapping features by percentage?

texts:
I have geometry records of type MultiLineStringZM and I am trying to INSERT new records that overlap some of the records currently in the table. If they were type MultiPolygon I could use ST_Area() as in the following

SELECT a.id
    FROM a, b
    WHERE ST_Area(ST_Intersection(a.geom, b.geom)) > 0.5 * ST_Area(a.geom);

to identify the polygons that overlap by a majority. However ST_Area is not a compatible function for MultiLineStringZM, and the above query does not return any results.

ST_Covers(), ST_Within() & ST_Intersects() will either not return any results or return too many depending.

How can I identify & delete all records overlapped by a 'majority' or some specified fraction for this geometry type?


texts: Shapefile containing several fields with the same variable at different times. How to group them into a single field?

texts:
I have a point shapefile that includes several (e.g 160) of fields with the mean temperature and wind speed [at each location] for different days. Ideally, I want to reduce the number of fields to only one for each variable, by grouping the corresponding data. For instance, what I was thinking to do, was to convert the corresponding values into strings and concatenate them into a long string, which I will store in a single field. What is the most efficient and elegant approach? Is it better to use a different data structure?

Please keep in mind that I have to use these data for plotting in matplotlib and also to be able to share the point shapefile and the temperature + wind speed data with others.


texts: Show loading icon with custom tool using ASP.NET and ArcGIS Server

texts:
I have a custom tool in vb.net for a web application, I'm using Asp.net, the class is in vb.net and the GIS service is in ArcGIS Server 9.3 sp1, the tool makes a selection of layers and takes a lot of time, so I want to show a loading icon on the map while it's getting the info from the server.
Any idea?
I already tried to modify the progress bar from the map but it works only within the map events, I want to show the progress bar from the custom class, from vb.net.


texts:
I am taking the context and taking it to a canvas which was asked here and the answer can be found in this JSfiddle example.

However, I am looking to load a WMS layer from geoserver and extract only this WMS layer to make a process and display something on canvas.
On Openlayers 2 it could be done with layerContainerDiv extracting the element that contains all the layers and taking the required one.
But it seems that in Openlayers 3 cannot be done. I have tried something similar to this example (take into account that this example gives a cross-origin error, but not in my server), but without success as all the layers are processed and cannot be distinguished.

Does someone know how can layers can be distinguished or separated to process them independently?

I found the Layer Spy Example which it seems that making use of precompose and postcompose could be a possibility, but  using the options of ol.Map I can find the 2 layers by map.getLayers().getLength(). And I can find the array of the WMS with map.getLayers().getArray()[0] but how to take the layer context?


texts: Where can I find geodata with postal codes of China with an open license (commercial use is allowed)?


texts: How to correctly import GeoJSON in PostGIS?

texts:
I'm trying to migrate GeoJSON data to PostGIS database and I thought that ogr2ogr is the tool:

ogr2ogr -f PostgreSQL PG:"host={} user={} port={} dbname={} password={}" test.geojson -nln ogr_test

However, the resulting PostGIS table contains geometry errors that were not present in GeoJSON file.

Using "Check geometry validity" from QGIS, shows that GeoJSON file is free of errors, but PostGIS layer that I converted with ogr2ogr contains multiple errors of this kind:

ring 3 of polygon 0 not in exterior ring
segment 2 of ring 11 of polygon 0 intersects segment 3 of ring 0 of polygon 0 at ...

with total 113 errors on 25 features.

Also issuing this query in pgAdmin:

SELECT id
FROM ogr_test
WHERE ST_IsValid(wkb_geometry) = 'f'

list 25 errors - the same number of feature errors as reported by QGIS command.

Then I converted my GeoJSON file to Shapefile. Testing the Shapefile for geometry errors, shows 9 errors on 2 features. Interestingly, the errors are on features that didn't produce errors while merged to PostGIS database.

Now I assume that ogr2ogr errors depend on destination format specifics, I'm interested as to why ogr2ogr produces these artifacts without informing the user, and what is the correct way to import GeoJSON to PostGIS?

I use GDAL 1.11.0

Actually it's again my fault - I used this script to check validity on GeoJSON files, as there were plenty of files to process. This script returned error only on one file which is not the test file used in this example. Loading GeoJSON in QGIS and checking validity, shows same errors.

So this seems my fault in arbitrary using multiple tools for validity check, and I'll vote this question to be close.

Sorry.


texts:
I am trying to make a tool for exporting pdfs from a directory: I found one online but it was very slow, I am trying to code it more efficiently. I am getting an AssertionError and I am not sure how to fix it.

import arcpy
import os

arcpy.env.workspace = r"C:\Users\jds\Desktop\MXD_Export_Test"

folderPath = r"C:\Users\jds\Desktop\MXD_Export_Test"

mxdList = arcpy.ListFiles("*.mxd")

for mxd in mxdList:
    saveName = mxd.split(".")[0] + ".pdf"
    savePath = os.path.join(folderPath, saveName)
    arcpy.mapping.ExportToPDF(mxd, savePath)

Error =

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.3\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in main.dict
  File "C:\Users\jds\Desktop\Scripts\Export_PDF_From_Directory.py", line 14, in
    arcpy.mapping.ExportToPDF(mxd, savePath)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\mapping.py", line 1148, in ExportToPDF
    assert isinstance(map_document, MapDocument)
AssertionError


texts: How to export a table as a geometry in Google Earth Engine Code Editor

texts:
I have imported my shapefile into GEE which is in my assets folder as a table. I extracted the shape from the table using this code:
table = table.geometry();
Map.centerObject(table);
Map.addLayer(table, {color: 'red'}, 'table');

This allows the shape to be visible but remains a layer and not geometry. I want to be able to have the shapefile in my assets readily available to import as geometry as the code I need to run requires the shape to be a geometry.


texts:
I want to extract road network data from OSM for a small city.

From where can I get the lane information of any road in OSM.

I tried exporting the data in the OSM XML format and there was no lane fields in that file.


texts:
I have a challenge viewing my file.

I get the following error (Failed: parse error at line 16032, column 25) when i try to open on Google.

I am unfamiliar with coding.


texts: how to improve poi visualization in tilemill

texts: I have imported all points of interest (POI) from database to TileMill. The distance between the POIs are a few meters. Thats is why its became ugly looking map. Is there any option that will aggregate POIs those are within 100 or 500 meter?


texts: Merge multilinestrings that are out of order

texts:
I am trying to merge a multilinestring object to create a single linestring object and it looks like the st_line_merge function is appropriate for this purpose. However, it only seems to work when the line segments are "in order".

Consider this first code block that makes a LINESTRING:

mls1 = st_multilinestring(list(rbind(c(0,0), c(0.5, 0.7), c(1,1)),
                              rbind(c(2,0), c(1,1)),
                              rbind(c(2,0),c(5,0))))
st_line_merge(mls1)

Compare with this block that does not create a LINESTRING:

mls2 = st_multilinestring(list(rbind(c(0,0), c(0.5, 0.7), c(1,1)),
                              rbind(c(2,0), c(1,1)),
                              rbind(c(0.5,0.7),c(-1,2))))
st_line_merge(mls2)

How can I make mls2 into a LINESTRING not a MULTILINESTRING?


texts: Remapping a raster map using GDAL in Python

texts:
I have a numpy array of data in Python that is a processed version (outliers set to zero, for example) of data from an Albers Conical United States (US) map:
PROJCS["Albers Conical Equal Area",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
I would like to take the data points and downsample them so that I can match them up with global weather forecasts, which are typically on an equirectangular projection with, say, pixels every 0.5 degrees.  So most of the globe will be empty -- there will just be some information in the US area.  I would like the information that is downsampled from the US Abers map to the equirectangular globe to be equal either to the average or the sum of the source pixels, depending on the data I'm downsampling.  The source pixels will, of course, not line up directly with the destination pixels.
I have figured out how to do this point-by-point, but it looks like gdal warp can do it for an entire map.  It seems like warp might even be able to do smart things with source pixels that span more than one destination pixel, but it's hard to tell.
Can gdal warp do this operation?  Can it do it from the Python interface?


texts: Selecting additional features in ArcMap 10

texts:
I'm having trouble selecting features in ArcMap 10. I have a shapefile of points, and I am in editing mode. I want to drag a box to select some points, then drag a new box and have these points added to my selection.

Holding down shift while creating the new box will invert the selection - so any previously unselected point are added to the selection, but any already selected points that are in the overlap are deselected. This makes selecting a large number of points difficult, as I must zoom in each time to ensure that there is no overlap.

How can I simply add points to my selection?


texts: QGIS prints proportional circles when value is 0 but not when NULL

texts:
Using QGIS3, I'm making maps with proportional circles based on an attribute value (integer).
When this attribute is NULL, no circle is printed, which is what I expect.
When this attribute is 0, a circle with size 1 is printed ; it should not.
I can't find any option to solve this.
In my sources, values are considered from 1 to 20 (max value). I tried with 0 as min value, it's the same.
There is an option to set size when value is NULL, but not 0.
I could set all 0 values as NULL, but i need the 0 to make operations (as sum).
Is there a way not to print circles on 0 values ?


texts: Determining features which are Near By Group using ArcGIS Desktop?

texts:
I have two polyline layers with corresponding attributes. Each layer has 100,000+ rows. I need to calculate the distance between the corresponding features from each layer.

E.g.

Layer 1:

FID | LineID
 1     1234
 2     5678
 3     9012

Layer 2:

FID | LineID
 1     9012
 2     1234
 3     5678

I need to find the distance between the rows from each layer that match by LineID.

ESRI had a blog about using model builder to create a Near By Group tool - https://blogs.esri.com/esri/arcgis/2010/09/16/nearbygroup/

The tool seems to stall out (probably because of the large number of rows). Does anyone have any other ideas on how I can accomplish this task?


texts: Privileges for Oracle users created using the ArcGIS Create Database User tool

texts:
I've created an enterprise geodatabase on Oracle. I can connect to it from ArcCatalog and from SQL Developer using the SDE user.

I've created a new user named GFPTest using the Create Database User tool, choosing the defaults for Role and Tablespace Name.

I can connect to the database in ArcCatalog using the GFPTest user. However, when I connect in SQL Developer using the GFPTest tool, I can't see any tables.

What have I missed in creating the database user that prevents me from viewing the tables in  SQL Developer?


texts: nDSM tree detection using scipy watershed

texts:
I am trying to generate segments with single trees in a normalized surface model (nDSM). I tried using scipy.ndimage to achieve this, but the segmentation results are not precise. I would like a better segmentation. How could I achive this?
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
import rasterio

from scipy import ndimage as ndi
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from skimage import filters

src = rasterio.open('nDSM.tif')
img = src.read(1)

img = np.where(img > 5, img, 0)
img = ndi.gaussian_filter(img, sigma=0.5)
distance = ndi.distance_transform_edt(img)
coords = peak_local_max(distance, footprint=np.ones((2, 2)), labels=img)

mask = np.zeros(distance.shape, dtype=bool)
mask[tuple(coords.T)] = True

markers, _ = ndi.label(mask)
labels = watershed(img, markers, mask=mask)


texts: Is there a tool to buffer features with more than one distance value in GRASS GIS?

texts: I´m looking for a tool like Multiple-Ring-Buffer in ArcGIS, but for GRASS GIS to be able to buffer with different distances but only having to start the process once.


texts: Creating Python Script Tool with Sortable Field Multi-Value Input?

texts:
I want to create a script tool in a regular ArcMap toolbox (not python tool box if possible) and have a multi-value field input parameter that allows you to sort the input fields, just like the "Attributes of Interest" input in the "Similarity Search" GP Tool in the built-in Spatial Stats Tools toolbox.  I can get a regular multi-value field input, but that only gives me all of the fields (that meet the filter requirements) in a pre-populated list, with check-boxes.  What my script tool needs however, is to get a list of field names from the user in a user specified order.

The "Similarity Search" GP tool has 2 inputs at the end of that tool that do exactly what I am wanting, and I've compared their settings with mine and I see nothing different, but they have a sortable multi-field input.  Additionally, if I copy that script tool into my own toolbox and delete the entire validation code section (code on the Validation tab of the script tool properties dialog), when I re-open that tool it still has the sortable multi-field input desired, so it does not appear to be validation code specific.  However, if I change any of the parameter display names in that script tool's properties, it suddenly reverts to the regular multi-field select box.

How can I get my script tool to replicate the Similarity Search GP tool's multi-field input?

I'd ideally like:

Instead of:


texts: How to add point features with coords taken from a .txt and show them on the map canvas?

texts:
Ok, so i am on Windows 7 and using QGIS 2.6 and i have a .txt file with some simple coords. I have extracted the coords to a list, no prob with that. But i looked at the PyQGIS cookbook on how to add features and the code is somewhat obsolete for example .addAttribute(0, 'hello') has changed to setAttribute. Anyway, i already have an empty point vector layer ready and 3 fields on it: id,x,y. I also need to fill them up as well, as i add the points to the map.
My code so far is:
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface

layer = iface.activeLayer()
points=[]
txt=open('C:/users/stelios/desktop/vasia/points.txt','r')

for line in txt.readlines():
    points.append(line.split())

print points
#so far so good, now i need to fix the part where i actually add the points and show them in the map

for x,y in points:
    feat = QgsFeature()
    feat.setAttribute(0, "1")
    feat.setGeometry(QgsGeometry.fromPoint(QgsPoint(float(x), float(y))))
    (res, outFeats) = layer.dataProvider().addFeatures([feat])

iface.mapCanvas().refresh()
iface.legendInterface().refreshLayerSymbology(layer)
txt.close()

When i try this i get a:
feat.setAttribute(0, "1")
KeyError: '0'
isnt the first value of setAttrubute, the index of the field i want to edit? So 0='id',1='x' and 2='y' or am i missing something? IF i remove the setAttribute line, then it works just fine and i see my points where they should be after i zoom to layer.


texts:
Is there a way to select what to import from planet-latest.osm.bz2 to database?

I want to do something like :

osm2pgsql ...  --place='city' --admin_level='2'

My goal is to use only cities and countries im my project ( maybe towns ).
I can download many databases with cities and countries shapes from the world, but I think OSM is more complete and if I decide to use other things is just import again.

I just try import all full planet data and then SELECT and INSERT in other database using SQL filters, but its too time-consuming ( and non-professional ).

Thanks.


texts: Getting SyntaxError when trying to loop through shapefiles in folder, calculate Euclidean distance and save files with appended name using ArcPy

texts:
I want to loop through shapefiles in a folder, calculate Euclidean distance and then save the files in an output folder with an appended name in arcpy. I want the name to be the filename with "KYU_" at the beginning and "_dist" at the end. Using @BERA code:
import arcpy, os

arcpy.env.workspace = "E:\\test_input" #There are two shapefiles in here
OutputFolder = "E:\\test_output"
arcpy.env.extent = "415637.469400001 2988625.7212 796637.469400001 3843125.7212"
arcpy.env.mask = "E:\\boundary\\boundary.tif"
for file in arcpy.ListFeatureClasses():
    newfile = os.path.join(OutputFolder, 'KYU_'+file.replace('.shp','')+'_dist.tif')
    print file
    print newfile
    arcpy.gp.EucDistance_sa(file, newfile, "550000", "250", Output_direction_raster, "PLANAR", "", Output_back_direction_raster)

Getting the error:
Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
TypeError: 'NoneType' object is not iterable


texts: ENVI5.0. Is it possible to use "Quickmap" option? How to create a simple map?

texts:
I used to work in ENVI 4.8, but now have to move on ENVI5.0
I found it was changed a lot comparably with the previous versions, obtained  many good new options, but some old options, in ENVI4.8, are disappeared or so significantly changed, that it’s not easy to find them quickly.

ENVI4.8 has a good simple option to create map from a scene rapidly – “Quickmap” (File > QuickMap > New QuickMap)
I can’t find this option in ENVI 5.
I can’t find any possibility to create map here (with legend, title, etc).
How is it possible to create a simple map in ENVI5?
Or, what option substitutes “Quickmap”?


texts: Order sql server spatial geometry distance by the nearest

texts:
I have a geometry column in a table, I want to get the N nearest neighbours; here is a specific example with the geography type

DECLARE @PTO Geography = Geography::STPointFromText('POINT(-57.432546 -25.365799)', 4326)
DECLARE @PDV3893 Geography = Geography::STPointFromText('POINT(-57.441357 -25.355435)', 4326)
DECLARE @PDV570 Geography = Geography::STPointFromText('POINT(-57.419036 -25.368309)', 4326)
SELECT @PTO.STDistance(@PDV3893) as n3893, @PTO.STDistance(@PDV570) as n570

This is the result of the query:

n3893               n570
1450,74600838172    1387,89423301988

Because the distance is in meters, you can see the nearest point is the n570 one

But in my query using the table (and therefore the geometry type) I get degrees, not meters

DECLARE @PTO geometry = geometry::STPointFromText('POINT(-57.432546 -25.365799)', 4326)
DECLARE @PDV3893 geometry = geometry::STPointFromText('POINT(-57.441357 -25.355435)', 4326)
DECLARE @PDV570 geometry = geometry::STPointFromText('POINT(-57.419036 -25.368309)', 4326)
SELECT @PTO.STDistance(@PDV3893) as n3893, @PTO.STDistance(@PDV570) as n570

Result:

n3893               n570
0,0136031693733486  0,0137411862661163

As you can see I can't perform a "order by" statement in a query with geometry because obviously we are working with degrees

How can I get the nearest point in geometry and order the results by distance?


texts: Distinct pattern symbols on the edges of overlapping polygons in QGIS

texts:
I have layers with different symbol fillings + - 0.

At the edge these symbols get cut off.

So I know this is possible for labels, but any way to only force full symbols?


texts: Check whether any feature classes in multiple gdb's has feature class representation

texts:
I am trying to find whether any feature classes in multiple gdb's has Feature class representation. If a feature class has represenation, "RuleID" & "override" fields will be in the attribute table. So i am just checking whether these two fields are present in the attribute tables.

This code is printing the results. But also returning few errors. Don't know whats the problem.

import arcpy

def FindField(fc,myField):
        fieldList = arcpy.ListFields(fc)
        for field in fieldList:
            if str.lower(str(field.name)) == str.lower(myField):
                print "    " + fc + " contains fieldname: " + myField

def FindField2(fc,myField2):
        fieldList = arcpy.ListFields(fc)
        for field in fieldList:
            if str.lower(str(field.name)) == str.lower(myField2):
                print "    " + fc + " contains fieldname: " + myField2

myField = "RuleID"
myField2 = "Override"
dir = r'D:\Test\gdb'
arcpy.env.workspace = dir

gdbList = arcpy.ListWorkspaces('*','FileGDB')

for gdb in gdbList:
    arcpy.env.workspace = gdb
    datasetList = arcpy.ListDatasets('*','Feature')
    fieldList = arcpy.ListFeatureClasses()
    for fc in fieldList:
        print arcpy.env.workspace,fc
        print "Searching root level Featureclasses..."
        print "  Searching " + fc
        FindField(fc,myField)
        FindField(fc,myField2)

    for dataset in datasetList:
        arcpy.env.workspace = dataset
        fieldList = arcpy.ListFeatureClasses()
        for fc in fieldList:
            print arcpy.env.workspace,fc
            print "  Searching Featureclass... " + fc
            FindField(fc,myField)
            FindField(fc,myField2)

D:\Test\gdb\RT_FactGeology_Qld_25K.gdb Rio_FactGeology_Qld_L
Searching root level Featureclasses...
Searching Rio_FactGeology_Qld_L
Rio_FactGeology_Qld_L contains fieldname: RuleID
Rio_FactGeology_Qld_L contains fieldname: Override

D:\Test\gdb\RT_Roads_Qld_25K.gdb Rio_Roads_Qld_L
Searching root level Featureclasses...
Searching Rio_FactGeology_Qld_L
Rio_Roads_Qld_L contains fieldname: RuleID
Rio_Roads_Qld_L contains fieldname: Override

Runtime error
Traceback (most recent call last):
File "<string>", line 31, in <module>
File "<string>", line 6, in FindField
UnicodeEncodeError: 'ascii' codec can't encode character u'\xf3' in position 7: ordinal not in range(128)


texts: Simple answer: Difference between a dataset and a map

texts: For a project I am involved in, we are creating a dataset that will show the extents of natural areas and provide an estimate of the ecological integrity of that area. We are calling it an "Ecological Integrity Dataset," but as far as I can tell, we're making a map layer that has this information. Can someone very simply tell me the difference between a map layer and a dataset?


texts: OGR2OGR EPSG:27700 Shapefile To EPSG:4326 MSSQL Table With Geography

texts:
I have been messing around with OGR2OGR for a while now trying to import a shapefile data into an MSSQL table. I am using the OGR2OGR executable that comes packaged with QGIS Wien (GDAL version 1.11.2).

I can get the shapefile to import when using the default Geometry data type however, when I try to output to a Geography data type (required for our GIS Reporting system), I get an output however, it doesn't seem to be able to display correctly in any of the systems that I use.

The command syntax that I am using is as follows:

C:\Program Files\QGIS Wien\bin\ogr2ogr -overwrite -f MSSQLSpatial "MSSQL:server=<My Database Server>;database=<My Database>;trusted_connection=yes" "<My Shapefile>" --config GDAL_DATA "C:\Program Files\QGIS Wien\share\gdal" -s_srs "EPSG:27700" -t_srs "EPSG:4326" -a_srs "EPSG:4326" -lco GEOM_TYPE=geography -lco SRID=4326

Here is what the output looks like in SSMS:

Any ideas?

I have a very limited tool set hence using OGR2OGR in this manner to start with.


texts: Integrating visual parameters into image file in order to export image

texts:
When I display the new image layer within the engine, it works like a charm. When I try to export the image file, however, it looses all of its properties that we are interested in for a simple CNN.

How can I include the properties of the visual parameters into the image file or into the export function (as there are no specific argument for them)?

The code:

//Landsat data
var geometry = ee.Geometry.Rectangle([-81.96, 32.64, -82.39, 32.94]);

var start = ee.Date('2014-01-01');
var finish = ee.Date('2014-02-28');

// I just use the dataset to get an overview over the available images for the given area.
var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate(start, finish)
                  .filterBounds(geometry);

var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

var image1 = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_017037_20140217');
//var image1 = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_017037_20140217')
//.visualize(bands(['B4', 'B3', 'B2']),  min(0),  max(3000), gamma(1.4));
//I have tried the two lines above but don't get it properly working. Any
//recommendations?

image1 = image1.toInt16();

Map.setCenter(-83.1406, 33.2426);

Map.addLayer(image1, visParams);

//What is the best way to carry the above specified information to the Export.image.toDrive() function?

//Export the image, specifying scale and region.
Export.image.toDrive({
  image: image1,
  description: 'Farms_Area_20140217',
  scale: 30,
  region: geometry
});

As this is just a short time project and as I have never programmed with JavaScript before, I am just looking for a fast and short answer!


texts: Ambiguous values on OSM to calculate optimal path: Which value should I use?

texts:
I'm using San Francisco OSM dataset, I want to calculate the optimal path using pgrouting(pgr_astar)
The problem that I don't know which column should I use?
For example:
I want to calculate the optimal path from Cesar Chavez Street to Tewksbury Avenue, which value of wich line should I use as source and target id?

Cesar Chavez Street:

Tewksbury Avenue:

For Cesar Chavez Street for example i have 115 street with the same name which one should I use and which value?


texts:
About a month ago we made the huge jump from ArcGIS 10.1 to 10.6.

Just recently we started noticing some odd behavior with graphics. Namely if you draw a graphic (say a polygon) on the map and then change the map's spatial reference the graphic disappears. Draw another polygon and then change the map's spatial reference and THAT graphic disappears. In short a graphic is only visible when its spatial reference matches that of the map.

I wondered if that was the way it always was so I went back and tested this on ArcGIS 10.1 and the oddity was not observed. Wondering if other users see this same issue in 10.6 or other versions between 10.1 and 10.

Fairly simple steps to reproduce:

Open ArcMap and add a country or continent feature class.
Set the map to wgs84.
Draw a graphic.
Change the map to Robinson World.
Draw a graphic.
Change map between wgs84 and Robinson.

Notice that the only graphic that is visible is the one generated in that spatial ref. The expected result is that all graphics should be projected on the fly. In ArcGIS 10.1 the results are as expected (ie all graphics get projected on the fly) but ArcGIS 10.6 graphics disappear when changing the map coordinate system.


texts: Question on how to use the editor tool and snapping in ArcGIS 10

texts:
I am new to ArcGIS 10 and am having trouble in using the new editing tool and snapping interface they have now. I cannot find the Task and Target options that allow me to choose Modify feature and define which layer I need to edit (and toggle between layers). Also, I could not find the Trace tool that used to be in ArcGIS 9.x versions that could be
useful for redrawing curvy boundaries. There seems to be a problem with Sanapping because
now there is a small toolbar that gives me only the options (End, Vertex,
Point, Edge, Tangent, etc) but without reference as to which layer I want
to snap to. Yet I can not snap properly to the vertex I want. Previously the cursor
would be snapped to the point I was hovering over within a reasonable proximity (defined by the tolerance). Now it seems that that sometimes does not happen even if I zoom in close enough. This of course results in problems since sometimes I got the snapping right, sometimes not.

Has anybody experienced similar problems and if so how did you resolve them? I would appreciate the information very much.
Thanks


texts: Geocode.farm addresses are not allowed to be displayed? Alternatives?

texts:
I just read the ToS of Geocode.farm and I am really curious about this paragraph :

  You may NOT display the information and data on either an internet or
  intranet site, regardless of access restrictions to such site unless
  reasonable precautions have been taken to hide and/or conceal the
  source data(coordinate pairs and/or address(es) returned from the
  Geocode.Farm API).

Do I understand this correctly - I can't use this service to e.g. convert GPS records to addresses that I show to the user?

This does not make any sense to me as this probably is one of the main use cases for reverse geocoding.

If this is true, does anybody know of an affordable and accurate alternative available for europe that allows storing and showing resulting addresses?


texts: Convert 32 bit image to 8 bit image in QGIS

texts:
I want to convert my images from 32 bits to 8 bits. I'm using otb to calculate texture and the raster output is on 8 bits, I'm going to use this information as a layer stack for my multispectral image at 32 bits!

How to convert from 8 bits to 32 bits, if it possible in QGIS?


texts: Creating feature dataset by pgAdmin (PostgreSQL) that can be useD in ArcSDE connection?

texts: How to create feature dataset in PostgreSQL (pgAdmin) that can by use in ArcSDE conection (ArcGIS)?


texts: Losing raster quality when using raster calculator (ArcGIS Desktop 10.5.1)

texts:
I'm currently working with DEM rasters that I hillshaded. The data I received contain nodata borders, probably the footprint of the originial mosaic dataset. As I received them separately, I can't work on the footprint to delete the nodata lines.
I'm trying to calculate the noData cells using the average value of the cells next to it, using raster calculator and this expression:

 Con(IsNull("DEM"), FocalStatistics("DEM", NbrRectangle(5,5, "CELL"), "MEAN"), "DEM")

The expression is working, but I'm losing quality of my hillshaded DEM.
Do you know how I can avoid this?

Before raster Calculator

After calculator on the down DEM


texts: How to fix a shapefile reprojected from sinusoidal projection to lat long?

texts:
I have been work in convert the MODIS Tiles Grid from sinusoidal projection to lat long projection using python ogr2ogr.

It doens't works well, the shapefile reprojected presents some lines that was not present in the original shape and it produces many mismatching.

Link to donwload the original MODIS Grid Tiles

Any suggestion to fix it?

# my code in command line
ogr2ogr -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" -t_srs EPSG:4326 -f "ESRI Shapefile" "C:/path/tiles_mod.shp" "C:/path/tiles.shp"

The shapefiles reprojected

The original shapefiles in sinusoidal projection:

Processed shapefile using ogr2ogr.


texts: Georeferencing Tiffs after manipulation in Photoshop?

texts:
I am exporting very high-res (7000x7000px) hillshades from ArcMap as tiff files. I am then doing some work in photoshop to the images.

How can I easily georeference these tiffs after working on them in PS?

Is there a way I can simply define the coordinate system and the bounds of the tiff without using ArcMap?

I am trying to avoid ArcMap because I seem to loose all me edits from PS and ArcMap just has a tendency to foul things up once I reimport. I have QGIS installed if that might be a better approach.


texts: Get list of editable fields while also adding 'SHAPE@' to the end

texts:
I need to get a list of field names that only includes editable fields and not the geometry field, but then add the geometry field back to the end of the list as SHAPE@ token.  I have code that runs successful, but I am horrible with List Comprehensions and I am sure there is something better out there.

update_fieldsED = [f for f in fields if f.editable]
update_fieldsNoShape = [fNS.name for fNS in update_fieldsED if fNS.type not in ["Geometry"]]
update_fieldsNoShape.append("SHAPE@")

Is there a more compact way to achieve this?


texts: WPS descriptions via OWS describeprocess

texts:
When following these instructions to obtain descriptions of various WPS descriptions it's specified to go to:

http://localhost:8080/geoserver/ows?service=WPS&version=1.0.0&request=DescribeProcess&identifier=JTS:buffer

On my geoserver instance though, that request yields the following response:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://mygeoserver:80/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="InvalidParameterValue" locator="service">
<ows:ExceptionText>No service: ( ows )</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

I do have wps installed and have successfully ran some wps requests. Is there something I need to do to enable ows?

Edit

I'm on version 2.10.0 and the request was edited to say localhost only to hide the real server name.
Also, http://localhost/geoserver/ows?service=wps&version=1.0.0&request=GetCapabilities does work

Here is a snip of capabilities:

Actually I just reran http://localhost:8080/geoserver/ows?service=WPS&version=1.0.0&request=DescribeProcess&identifier=JTS:buffer and it's working now. Not really sure what was going on yesterday. I'll keep an eye on it; what's the proper stack etiquette in this case? Do I just close the ticket?


texts: How to force deletion of corrupt ArcSDE Raster?

texts:
ArcSDE 9.3 / ArcGIS Server 10.0 / ArcCatalog Advanced 10.0

After an unsuccessful attempt at creating a mosiac, ArcCatalog refuses to allow deletion of the corrupt Raster file.

The Mosiac creation was corrupted because of lack of file space on the MS SQL Server. Now I cannot open the file or delete it.

Is there any way I can force delete the Raster?


texts:
I use a android app called OSMTracker. This app records the GPS track and can take pictures and notes. Using JOSM software, the photos are shown exactly at their coordinates (location) but i cannot use these photos in other GIS softwares like QGIS, it says no geotagged photos in the location.

How can i find out whether a photo is geotagged or not? I know that the information is written in some EXIF header, but i don't have any other files. I only have the GPS track and the photos.

If the photos are not geotagged, can i geotag them using the GPS track?


texts: Script to Compute the river slope with correct DEM elevations?

texts:
I know this topic has been discussed before in this forum , I read all the suggestions but I think I still have the same issue.

Working with QGIS, I'm supposed to compute the slope of river channels and then be able to classify them in classes. For that, I'm working with a DEM 90m resolution and a river network derived from the openstreetmap database.

As my DEM gave me wrong values for some slopes because it sometimes don't consider my rivers as lowest points of my valleys, I would like to run a script that enable me compute my slope using obviously the right elevations of the DEM.I've joined a picture of the network I'm working on and here are the manipulations I've done so far after many days of vain efforts.

I decided to split my streams with the extension "Network" into segments with different lenghts and using the tool "connect" I've got nodes with IDs that increases from upstream->downstream .
I assigned to each of these nodes its corresponding elevation from the DEM 90m

In the picture black numbers corresponds to the ID node which are classified and the red numbers correpond to the elevations!

As you can see I've circled  two nodes where the algorithm I guess don't recognise the bifuraction of the river! I wonder if there's a solution to this issue too!

I've got sometimes wrong elevations at some points ( e.g the elevation sometimes increases going downstream!!) I guess it's normal because this odd point sometimes is the value of a pixel located on a tree for example in the bank river! so the correct elevation is just in the adjacent pixel !

Being not well skilled in coding, how can I go from a point in the river to the  upstream and whenever there's an odd value of elevation, the script skipped it and goes to the right one upstream and then calculate the distance between the points so as to find the slope?


texts: Getting started with development for ArcGIS Server

texts:
I'm currently working on upgrading a C# web service I wrote a couple of years ago built against ArcGIS Server to ArcGIS Server 10. The application in question does some basic spatial queries and exports several PDFs, so no addtional extensions are required however it cannot be hosted directly on the ArcGIS Server.

My problem is that I'm unsure of exactly which ESRI components I'll need and where they need to be installed.

Is what I've drawn up below reasonable?

Developer PC

Visual Studio 2010
ArcEngine 10 Runtime
ArcObjects 10 SDK

Application Server

ArcEngine 10 Runtime

ArcGIS Server

ArcGIS Server 10
ArcGIS GIS Services


texts:
in google map,

how to show Douglas county MAP in regions manner. (east, west, north, south).

These regions should cover with border and show the picker when click the balloon.

Here is my KML file, which shows the Douglas county, now i want to split under region.
(East,West,North,South, NE,NW,SE,SW)
my kml file

this kml output display here.

i posted this post in stackoverflow, but it was closed due to ambiguous. can you please some share your ideas.


texts: In R The Volcano dataset has elevation values, but is not a DEM file. How can I convert my DEM files (from earthexplorer.usgs.gov) so that the I can access all the elevation values in a matrix (what is returned by volcano)? I am trying to use r2stl() with the data.


texts: How to SUM field and then insert into NULL column?

texts:
I am only in the beginning stages of learning Python.

I have a column, 'QUANTITY_SOLID', that I want to take the sum of all of its rows. I understand I can do this with 'SUMMARY_STATISTICS', but i'm not trying to create another entire table nor do anything other than simply SUM up that column. Once I have it summed up, I want to insert the sum into a new NULL field that I have created - I am attempting to do this via the 'CALCULATE_FIELD' tool.


texts: Calculating NDVI values using Sentinel 3 OLCI

texts:
I am trying to compute NDVI values for certain polygons using SEN3 dataset, but keeps getting negative values
It seems I am getting wrong values for NIR and RED bands when trying to transform SENTINEL 3 OLCI values from radiance to surface reflectance.
I am assuming that is the problem, maybe not?

var sigma2VisParmSEN3 = {min: 0, max: 1, palette: 'red,orange,yellow,green'};

var SEN3 = ee.ImageCollection('COPERNICUS/S3/OLCI')
                  .filterDate('2019-01-01', '2019-02-01')

// Select bands for visualization and apply band-specific scale factors.
var SEN_corr = SEN3.select(['Oa17_radiance', 'Oa08_radiance'])
              .median()
              // Convert to radiance units.
              .multiply(ee.Image([0.00493004, 0.00876539]));

print(SEN_corr);
//calculate NDVI values per region - ***please change region to a shape you can use on your device***
var ndviSEN3 = SEN_corr
      .clipToCollection(region)
      .normalizedDifference(['Oa17_radiance', 'Oa08_radiance'])
      .rename('ndviSEN3');

Map.addLayer(ndviSEN3, sigma2VisParmSEN3, 'ndvi2SigmaSEN3');
var ndviFeatures = ee.FeatureCollection(ee.List(ndviSEN3
  .reduceRegion({
    reducer: ee.Reducer.toList(),
    geometry: region,
    scale: 10
  })
  .get('ndviSEN3'))
  .map(function (ndviSEN3) {
    return ee.Feature(null, {ndviSEN3: ndviSEN3})
  })
)
print(ndviFeatures.aggregate_array('ndviSEN3'))

I get these values for ndviFeatures:

(-0.06449873745441437, -0.3388196527957916.......)

Any suggestions?

Thank you!


texts:
I'm following these instructions from the FCC: https://transition.fcc.gov/form477/Geo/finding_census_blocks_in_qgis.pdf

Everything works just fine up until step #6, where the Intersect tool is used. I do everything just like the instructions, but when the Intersect tool finishes, it's just added an empty layer. I've hidden all the other layers, so it's not just getting buried underneath them or anything.

Any ideas as to what might be going wrong?