texts: Failed Create Space Time Cube from Multidimensional Raster tool

texts:
I am creating a space-time cube from global forest watch forest loss data. I have created a multidimensional raster of years from 2001-2020 of forest loss cells. I run the Create Space Time Cube from Multidimensional Raster tool in ArcGIS Pro 2.7.3. I ran from the geoprocessing pane, I used the default settings with view to adjust after the first successful run. I used an multidimensional raster in .crf format meeting the conditions required including 10 time steps. The tool fails and I got this message:
Traceback (most recent call last):
  File "<string>", line 4594, in execute
  File "c:\program files\arcgis\pro\Resources\ArcToolbox\Scripts\SSCube.py", line 92, in __init__
    self.__initialize(cubeObj)
  File "c:\program files\arcgis\pro\Resources\ArcToolbox\Scripts\SSCube.py", line 168, in __initialize
    self.__createDimensions(cubeObj)
  File "c:\program files\arcgis\pro\Resources\ArcToolbox\Scripts\SSCube.py", line 540, in __createDimensions
    self.__createLocationVariable(lat, 'latitude', latValues)
  File "c:\program files\arcgis\pro\Resources\ArcToolbox\Scripts\SSCube.py", line 626, in __createLocationVariable
    var[:] = value
  File "netCDF4\_netCDF4.pyx", line 4992, in netCDF4._netCDF4.Variable.__setitem__
  File "netCDF4\_netCDF4.pyx", line 5271, in netCDF4._netCDF4.Variable._put
  File "netCDF4\_netCDF4.pyx", line 1928, in netCDF4._netCDF4._ensure_nc_success
RuntimeError: NetCDF: HDF error
 Failed to execute (CreateSpaceTimeCubeMDRasterLayer).

What is causing the tool to fail and how can I fix it?


texts: Extracting a polygon from slanted Sentinel satellite image using R

texts:
I have a Sentinel image, the coordinates of the boundary are given. I want to extract/subset a polygon from it. I am using R and don't have any prior experience in GIS.

I have downloaded some Sentinel 1 images from Copernicus https://scihub.copernicus.eu/dhus/#/home. The associated KML file has the coordinates of the image corners.

<coordinates>-118.362823,34.353683 -121.078079,34.752659 -120.786476,36.256836 -118.020294,35.860096</coordinates>

If I plot the coordinates on a map using DaftLogic Area Calculator it looks like the screenshot below:

The shape drawn above isn't parallel/perpendicular to the latitude and longitude lines, its inclined because of the satellite orbit. Also, two of the sides have a difference of ~1 km. So the shape is not in rectangular format.

I am opening the satellite image using raster package in R. The extent(georeferencing) of the image/raster doesn't have the actual geo coordinates. But, it represents the dimensions of this image.

library(raster)

s <- brick('s1a-iw-grd-vv-20161030t230720-20161030t230745-013726-016066-001.tiff')

s
#class       : RasterBrick
#dimensions  : 16716, 25485, 426007260, 1  (nrow, ncol, ncell, nlayers)
#resolution  : 1, 1  (x, y)
#extent      : 0, 25485, 0, 16716  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#data source : s1a-iw-grd-vv-20161030t230720-20161030t230745-013726-016066-001.tiff
#names       : s1a.iw.grd.vv.20161030t230720.20161030t230745.013726.016066.001
#min values  : 0
#max values  : 65535

extent(s)
# class       : Extent
# xmin        : 0
# xmax        : 25485
# ymin        : 0
# ymax        : 16716

There is a function to assign the extent, setExtent(). This would have worked if the shape/polygon was rectangular and not inclined.

What I need to do after the correct georeferencing is done, is to get a subset of the raster based on the polygon using extract(). Since I am not able to set the extent properly, I am not able to do so.

Why R?
I can use any other language as well. I am using a language because I have to repeat the process with multiple images.

One more problem to deal with is, some of the images are mirrored image, this adds to the problem. Note for mirrored image https://sentinel.esa.int/web/sentinel/faq

Is this the right approach? Should I work with other package or language? If somebody has already done it, please help.


texts: How can I preserve image quality when combining aerials of different resolutions?

texts: I have two different rasters that were flown at different times during the same month. One has a resolution of 4-inches, the other has a resolution of 6-inches. Using ArcGIS tools, how can I combine these rasters while preserving as much image quality as possible? I intend to use the combined output as a basemap.


texts:
I have a road shapefile of a selected area of 1km x 1km created in QGIS. I have also been able to create a grid shapefile of 1m x 1m intervals over my area as below: -

Is there a possibility to get each grid cell's data (i.e, "does this cell contain a road or it is not a road?" kind of data taken from the road shapefile and the area shapefile that is in its background) into a CSV file?


texts: ApplySymbologyFromLayer_management sets polygon outlines to default

texts:
I'm struggling with an issue which I think is a bug in arcpy.  I am using AddSymbologyFromLayer_management to apply the symbology from a reference layerfile (yes, an actual layerFILE, not a layer within the .mxd) to a newly created polygon shapefile which is displayed as a layer.

This works fine, even when I change what the field the symbology shall be based on. However, when I also change the class breaks (I use GRADUATED_COLORS here), the outlines, which are set to "No color" in the reference layerfile, are set back to default, which means some kind of grey with 0.4pt width.

This behaviour only occurs when I change the class breaks via code, if I keep them as they are, everything is fine.

This is what I basically do within the code:

# get reference layer
ref_lyr_file = os.path.join(workspace, 'styles', 'graduated_colors_reference_style.lyr')
# get newly created layer
new_lyr = arcpy.mapping.ListLayers(mxd, lyr, mapDF)[0]
# import symbology from reference layer-file
arcpy.ApplySymbologyFromLayer_management(new_lyr, ref_lyr_file)
# set value field
new_lyr.symbology.valueField = draw_param_dict[draw_parameter]
# reclassify values
new_lyr.symbology.reclassify()
breaks = new_lyr.symbology.classBreakValues
# round values nicely as new class breaks
lenMin = len(str(int(breaks[1])))
new_breaks = [int((round(b / (10 ** lenMin), 2)) * (10 ** lenMin)) for b in breaks]
# update map document
new_lyr.symbology.classBreakValues = new_breaks
arcpy.RefreshTOC()
arcpy.RefreshActiveView()

How I create the new breaks is not important, so you can skip this section. The thing is: the polygon outlines are set back to the default color ONLY AFTER applying new class break values.

Does anyone know how to get around this?


texts: QGIS - How do you keep a map fixed to a template?

texts:
I am creating templates within print composer for people within my company to use. I would like them to include a map of the county in the top left corner that can then be used as an overview for whatever map they have created. How can I set it up so that they can just open the template with the map included without them having to open that map? When I add items from template and then save the whole thing as a template and open it later from a new project it opens up blank. The image I add in the top right corner also doesn't remain. I have locked them when I have saved them.

What am I doing wrong?


texts: Classifying raster that takes neighbours into account?

texts:
I am working on classifying forest types using Sentinel-2 satellite imagery (10-m resolution; band-2, band-3, band-4 & band-8). I have ground-truthed field data as a training set. I want to set a parameter that constrains classification so that nearby pixels must be classified into the same class or that there must be a minimum of 10 neighbouring pixels within the same class. Ideally, I want the classification to take into account the heterogeneity of an area and not just classify each pixel based on the average in the training set. I would prefer a solution in R but I could also use QGIS or ArcGIS.

Details:

I have been using this helpful tutorial for the classification in R:
http://ceholden.github.io/open-geo-tutorial/R/chapter_5_classification.html.

While I have successfully run the functions and algorithms to produce a classification using the randomForest package and the raster package, I am not satisfied with it because the result is lots of isolated pixels of different classes. (For example, five of the six classes are adjacent neighbours with no dominant class to use in post-processing to cluster together). Visually, I noticed that some forest types are distinguished by being heterogeneous (pixels of very different colours side-by-side) and other forest types are homogeneous (pixels of similar colours side-by-side). I want the classification algorithm to pick it on the pattern within the training set (polygons of circles 60m in diameter) and/or to take into account neighbours in the classification so that nearby pixels have the same value. I also have polygons of delineated forest types from photo-interpretation, which I could use to constrain to have only single value per polygon (the values of the photo interpretation are not accurate but the limits are good enough).

I've taken screenshots of the RGB true-colour image of three (out of six) of forest types and their respective classifications. I've chosen two for which it works well and I could use post-classification majority filter to smooth groups. And one for which there is a high heterogeneity and the classification algorithm does not work.

What function can I use in R to classify raster that takes neighbours into account?


texts: PostGIS Connection In QGIS. Vector Query Results of GeometryCollection Load Layer as Raster

texts:
https://postgis.net/docs/ST_VoronoiPolygons.html

I am trying to load this query result using PostGIS 2.5 connection in QGIS 3.4, QGIS 2.18.

  SELECT
      ST_VoronoiPolygons(geom) As geom
  FROM (SELECT 'MULTIPOINT (50 30, 60 30, 100 100,10 150, 110 120)'::geometry As geom ) As g;

Getting the error of the vector layer GeometryCollection result as Raster Layer in QGIS.

Results =>

I have tried these PostGIS functions to verify results,

https://postgis.net/docs/ST_IsValid.html
https://postgis.net/docs/ST_MakeValid.html
https://postgis.net/docs/ST_IsValidReason.html

I have also applied ST_AsGeoJSON and ST_AsText, results is GeometryCollection, and valid.

I have imported the GeoJSON in https://mapshaper.org/, it's okay.

After extracting/dumping the Features from GeometryCollection with following functions it's working fine.

ST_CollectionExtract with Polygon selection.
ST_Dump.

Why QGIS Load GeometryCollection vector query result as a raster layer?
Is there any method available under QGIS to load GeometryCollection?


texts: Iterating over all features in featureclass to intersect them individually with polygon using ArcPy?

texts:
I've got a point featureclass, and I'm trying to individually intersect each feature with a polygon featureclass to get some information about which polygon each feature is in.  I did one intersect with an entire point featureclass against a polygon featureclass, which was straightforward, but I'm not grasping the model for doing the individual intersections.  First attempt below:

# -*- coding: utf-8 *-*
# Name: IntersectingPolygon.py
# Description: Find intersecting polygon on a feature by feature basis

# Import arcpy and set path to data
import arcpy
import csv
arcpy.env.workspace = "C:/dev/FGDBs/test.gdb"
intersecting_polygon_layer = 'pg_lyr'
feature_layer = 'feature_lyr'
test_feature = 'some_feature'

# Make a layer of the polygons
arcpy.MakeFeatureLayer_management('polygons', intersecting_polygon_layer)

feature_rows = arcpy.SearchCursor(test_feature)

count_features = 0
for feature_row in feature_rows:
    # This line fails, I think because you can't make a layer this way
    arcpy.MakeFeatureLayer_management(feature_row, feature_layer)
    arcpy.SelectLayerByLocation_management(intersecting_polygon_layer, 'intersect', feature_row)
    polygon_rows = arcpy.SearchCursor(intersecting_polygon_layer)
    print 'Feature Globalid:',feature_row.getValue('globalid')
    for polygon_row in polygon_rows:
        print '\tPolygon objectid:',polygon_row.getValue('objectid')

    if count_features > 5:
        break
    count_features = count_features + 1


texts: Calculate Galileo satellite position from Ephemeris data

texts:
I'm trying to calculate Galileo stallite's position from ephemerides data.

I'm reading from https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf at 3.3.1.

The first step is:

  Compute the time Tk from the ephemerides reference epoch toe (t and
  Toe are expressed in seconds in the GPS week):
  Tk = t - Toe

I have the value of Toe as it's retrieved from ephemeris data, but I'm not sure where the value of t comes from.

I have a hunch that t is retrieved fom the clock parameters, but those arn't in the ephemerid parameters. So if I retrieve the clock parameters separate from the ephemeris parameters then they will have the be retrieved at different times and create calculation errors errors.

I've also read https://www.navcen.uscg.gov/pubs/gps/sigspec/gpssps1.pdf at 2.5.7 which explains t a little bit:

But I don't understand what it means by the "effective satellite PRN code phase time at the message transmission time" and "PRN code phase time offset".

EDIT:

I haven't found a direct solution but I've found a different path to take so I'm gonna post it here in case someone else finds this useful. There is a Volume 2 of the first PDF I posted that contains more practical examples of the math. The second volume contains a CD that has a program called gLab Tool Suit with all the examples and other nice programs. One of the examples is "sub_orbit.f" which is a function that takes in the Ephemeris parameters and Epoch and returns the satellites position. I've created a pastebin for the example code here: https://pastebin.com/9b38LB1R (It's in Fortran) but you can download both the PDF and CD for yourself here: https://gssc.esa.int/navipedia/index.php/GNSS:Tools


texts: Exporting crown boundaries from tree segmentation in R?

texts:
I am using lidR package to segment my trees using lastrees directly from the points clouds.

The function works fine, however I would like to export to vector file the canopy boundaries from tree_segmentation that I did directly from my point cloud.

Is it possible to extract canopy boundaries from segmentation lastrees in order to analyze my results I got?

las = lastrees(lasn,li2012(R = 3, speed_up = 5))
plot(las, color = "treeID", colorPalette = col2)


texts: Summarizing values of polygons into other polygons using PostGIS?

texts:
I have two tables. One for school districts, one for schools with the total number of each student in the school.

I created a float field in the districts table for student population. I want to populate this field with the sum of students in each school within each district.

Essentially, I need to select by location the schools in each district and assign the summary of their student population to the district student population field.

I feel like this is a simple task with ArcMap, but I can't afford their license.


texts: Enabling model tool even if field already exists using ArcMap ModelBuilder?

texts:
I made a quick little tool in ModelBuilder to add specific fields based on a precondition (checked boxes). The issue I'm having is that if any of those fields already exist in the feature I'm trying to add fields too, the tool won't run. Sometimes some of the fields already exist in a feature, and I'm wanting to just add the missing ones. I could delete the existing fields, but then I would lose the data in those fields.

I've found info on a If-Field-Exists logical tool, but I can't seem to find it within Model Builder.

Is this a tool found only in ArcGIS Pro?


texts:
My problem is that I'm trying to extract contours and elevations from a DEM map. Everything else looks fine, but I think that the "ELEV" values I'm getting are not the real elevations, maybe I'm missing something. The steps that I follow are:

Import WMS DEM map to my project
Run the map to raster algorithm - [Settings: Map units per pixel: 1 (I'm doing tests with small areas) / make the background transparent / unique layer to render: the WMS DEM layer)
With the output I get from the first step, I run the "Contour" algorithm with the following settings: Band(Red) / Interval between contour: 1 (here I always use the same value as the "map units per pixel" value, correct me if I'm doing it wrong)
And that's it, the last thing I do is add the single labels that show the "ELEV" values of the contours.


texts: Output from selected polygons and source layer in a table

texts: How to output to a table a specified attribute from selected polygons and the source layer's line feature name or ID that intersects each selected polygon?  Basically trying to get a table that lists each well and the expiration dates for each leased parcel a well crosses.


texts: How do you use the setErrorMessage(message:string) within the ToolValidator Class?

texts:
The Programming ToolValidor Class ArcGIS Help menu states that: the following method:
setErrorMessage(message:string)
"Marks the parameter as having an error (a red X) with the supplied message. Tools do not execute if any of the parameters have an error."
Does anyone have an example of the code working to make the red x appear?


texts: Is there a Leaflet for Dummies or Leaflet 101

texts: I am very interested in learning leaflet, but I am not a programer nor a cartographer. Is there any LEAFLET for DUMMIES course that I can take?


texts: earthquake effects on total station surveys

texts: What I have been wondering is whether an earthquake can affect a closed system survey.  We are using Leica TCR 307s and TCR 407s and I was wondering if an earthquake can throw off your fixed points if staying within a 500m grid?  The fault line does not run through the site.


texts: How can I add output vector to a list to loop a geoprocessing on it?

texts:
I have about 18000 circular buffers that I want to transform in some kind of a wind rose. To do this, I transform the buffer polygons to lines, create points along the lines on the wind directions that I'm interested in (N, NW, W, SW, S, SE, E, NE), create voronoi polygons from theses points then intersect them with the original buffer polygons.

The problem is that some polygons are superposed, which breaks the voronoi processing.

To resolve it, I split the buffer layer into as many layers as I have polygons (so about 18000), and I want to loop my processing on each layer and then merge them into one.

To do this, I wrote this code :

#creation_rose_vents=name
#buffers=vector
#buffer_vent=output vector

import os, processing

path = "H:\\donnees_geo\\buffers_individuels"

shpfiles = [os.path.join(d, x)
            for d, dirs, files in os.walk(path)
            for x in files if x.endswith(".shp")]

list_buffers_vent = []

for buffers in shpfiles:
    outputs_QGISPOLYGONSTOLINES_1=processing.runalg('qgis:polygonstolines', buffers,None)
    outputs_QGISCREATEPOINTSALONGLINES_1=processing.runalg('qgis:createpointsalonglines', outputs_QGISPOLYGONSTOLINES_1['OUTPUT'],785.0,0.0,5495.0,None)
    outputs_QGISVORONOIPOLYGONS_1=processing.runalg('qgis:voronoipolygons', outputs_QGISCREATEPOINTSALONGLINES_1['output'],0.0,None)
    outputs_QGISINTERSECTION_1=processing.runalg('qgis:clip', outputs_QGISVORONOIPOLYGONS_1['OUTPUT'],buffers,buffer_vent)
    list_buffer_vent.append(buffer_vent)

So far it is just supposed to add the result of my processing to a new list (for which I will try to merge all objects, but this is another problem), but there is an error saying : "buffer_vent is not defined" and I don't understand what is wrong.

Does anyone see what I'm doing wrong?

Bonus question: Is it possible to merge all objects from a list to one single shapefile or should I write my output to shapefiles first and then merge them?


texts: Adding feature class within file geodatabase in QGIS with Python?

texts:
I've been using ArcPy for quite a while, but I'm now trying to learn how to create geospatial workflows using PyQGIS. I'm starting with simple things, like adding a layer to the map. I'm trying with the following code:

import os

gdb_URL = r'...'
fc_path = os.path.join(gdb_URL, 'feature_class')
layer = iface.addVectorLayer(fc_path, 'layer_name','ogr')

That throws a message saying that 'the layer is not valid'.

=> If I use the same piece of code to add a shapefile instead of a feature class within a file geodatabase, the layer gets added without problems.

=> If I try to add the feature class within a file geodatabase using the "Open Data Source Manager" button in QGIS, the feature class gets added without problems.

Any ideas why that piece of code is not valid to add feature classes?

Maybe 'ogr' is wrong and I should be using a different provider name.


texts: Creating service area using both streets and buslines in same ArcGIS network?

texts:
I want to create a network analysis more difficult. Therefore, I want to calculate a Service Area around hospitals next to bus stops. But: The network does not only consist of roads and pedestrians (time) but also of bus lines and bus stops. In best case the result should be a augmented service area, which includes bus drives to the hospital in a certain time! I nearly have the solution, but here is the problem. Trying to solve the analyse I realized that several bus stops are not within the calculation. I have checked the topology of the stops and lines and there are no errors (as ArcGIS says). Also the attribute table seems to be correct. I tried to planarize the bus lines but I can't use the tool, as it is not possible to click on it (grey color).

Maybe you have a idea to help me?

Problem picture:

(Source)


texts: Statistics of each image in an ImageCollection in Google Earth Engine

texts:
I have an ImageCollection that I made calculation with so that it now contains images that show a ratio of one band within a given Polygon. I would like to get the mean and standard deviation for each image of the collection within the polygon, so I have a value for each image afterwards. I keep on receiving the error of ReduceRegion not working for all the codes I am trying. Maybe you can help me to understand this better and solve the problem.
I tried this code for example:
I edited code but unfortunately still the same error. I read it could be because of null values but I do have a function sorting out null values so that shouldn't be the issue
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

function stats (image) {
  var stats1 = image.reduceRegion({
    reducer: reducers,
    geometry: MODISofstation,
    scale: 30,
    bestEffort: true
  });
  return image.set(stats1);

var statistics = ImageCollection.map(stats);
print (statistics);


texts: Add a OpenStreetMaps road data layer on Google Maps API v3

texts:
I've successfully embedded a Google Maps API v3 map on my site, and added OpenStreetMaps (OSM) as the base layer. Is it possible to use Google's tiles as the base layer, and then overlay street data on top of that, rather than the entire OSM layer?

<script type="text/javascript" charset="utf-8">
    $(function(){

        google.maps.MapTypeId.push("OSM");

        var osmMapType = new google.maps.ImageMapType({
            getTileUrl: function(coord, zoom) {
                return "http://tile.openstreetmap.org/" + zoom + "/" + coord.x + "/" + coord.y + ".png";
            },
            tileSize: new google.maps.Size(256, 256),
            isPng: true,
            alt: "OpenStreetMap layer",
            name: "OpenStreetMap"
        });

        $("#map-canvas").gmap3({
            action: 'init',
                options:{
                  center: [22.49156846196823, 89.75802349999992],
                  zoom:2,
                  mapTypeId: google.maps.MapTypeId.ROADMAP,//google.maps.MapTypeId.OSM,
                  mapTypeControl: true,
                  mapTypeControlOptions: {
                    style: google.maps.MapTypeControlStyle.DEFAULT
                  },
                  navigationControl: true,
                  scrollwheel: true,
                  streetViewControl: true
                }
        });
    });
</script>


texts: Vector layer won't move as window resizes in Chrome

texts: I have a GeoJSON layer that I am treating as a javascript object in my openlayers map. When I use IE or Firefox, everything works fine. However, in Chrome, as I resize the window, the GeoJSON layer is not moving with the background wms.It stays static and becomes displaced from the wms. Even when I click around the map, it stays in its incorrect placement. Any ideas why Chrome is doing this? I guess I could redraw the layer everytime the window is changed. But that seems too drastic for this problem with this browser. Thanks for any help!


texts: ExecuteError: ERROR 999999: Error executing function. Failed to execute (CreateFileGDB) ArcGIS Server 10.5

texts:
I am trying to convert geopackage file to geodatabase using my python script. It fails with 999999 error. This gdb is further used for the arcgis rest map service. It is only working when the service is stopped. Also, the schema lock is unchecked in the arcgis rest map directory. There is no lock on the gdb folder and no other Esri product is running on that server. I am not sure if it is related to the permissions on the server. Any suggestion is welcome.  Using Python 2.7.12 installed with ArcGIS Server 10.5
import sys
import os
import arcpy
#Set workspace
arcpy.env.workspace = r"C:\myFolder\gdb"
# To allow overwriting the outputs change the overwrite option to true.
arcpy.env.overwriteOutput  = True

gpkg_file = r"C:\myFolder\gpkg\geopk.gpkg\filename.file"
# Geodatabase Name
gdb_name = "myGdb.gdb"
out_path= r"C:\myFolder\gdb"

# Process: Create File Geodatabase

arcpy.CreateFileGDB_management(out_path, gdb_name, out_version="CURRENT")
print("geodatabase " +gdb_name+ " is created")

# Process: Feature Class to Feature Class

out_path = r"C:\myFolder\gdb\myGdb.gdb"
#geodatabase file name
gdb_file_name ="indexarea"

#arcpy.FeatureClassToFeatureClass_conversion(in_features, out_path, out_name)
arcpy.FeatureClassToFeatureClass_conversion(gpkg_file, out_path, gdb_file_name)
print("***geopacakage file successfully converted into geodatabase***")


texts:
I have write the following code for showing coordinates of map but I getting this error:

Uncaught TypeError: map.getLonLatFromViewPortPx is not a function

Here is my code :
map.on('click',function(e){
    map.forEachFeatureAtPixel(e.pixel,function(feature , layer){
    var LonLat = map.getLonLatFromViewPortPx(e.xy);
        LonLat.transform( map.projection,map.displayProjection);
        Longitude = LonLat.lon;
        Latitude  = LonLat.lat;
        console.log(Longitude);
        })


texts:
There's also a degree, minutes, and seconds format, and one with a decimal number, but what is the name of this format?

What is the name of the GPS format of this lat/long: 35-39-24.720N / 110-26-06.720W

I've seen it in a few different places, but Google Maps doesn't seem to understand it.


texts: My boundaries shape file has the following attributes; CODE, PROVINCE, AREA, SHAPE_LENG and SHAPE_AREA (note that there are no geo- coordinates). My data on the csv file has the following fields; PROVINCE, MALE, FEMALE. What l need is to have proportional symbol like a circle representing the males and females for each of the provinces on the map. How can l accomplish this???


texts: Using field calculator on shapefile converted from DWG in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.3.

A third party convert complex autocad drawing (dwg) to shapefile for us. I got from them a zip folder with 1 to 16 Shapefile (depending of the original dwg data), all shapefile come with 4 file (.dbf, .prj, .shp and .shx) and I can visualize the data without problem, the attribute table seem fine when I open it.

The problem arise when I try to update a field using either the field calculator or manually.
The below picture illustrate the problem, first is the original table, second is after using the field calculator and third the calculator.

(for test purpose the calculate was only writing "A" in the ALT column (column in text format, length 255))

As you can see not only spaces are inserted before the A (as if I have "   A" as input) but the POINT_X field got truncated and the POINT_Y and POINT_Z column are filled with what seem to be random value.
The same thing happen if a manually edit a value, as soon as I save edits all value are modified.

the insertion of space before text column follow a pattern : no space added on first column, one space added on second column, two space on third and so on...

This kind of behavior happens on all field (text field get added space and numeric field are modified) on all shape from this provider. The only workaround is to export the shapefile (exported shapefile come with .cpg, .dbf, .prj, .sbn, .sbx, .shp and .shx and behave as expected).

I will ,of course, ask the provider of the shapefile to look into the problem and try to fix it, but I'm curious what it could cause this behavior and how to fix it on my side.

adding a .cpg file from a working shapefile doesn't change anything but exporting either to .shp or to gdb fix the problem.

This problem seem to be a specific Arcgis problem as the problem didn't occur on QGIS (can edit without problem) and opening a messed-up shp in QGIS show the text column without the extra space before the text but the numeric column are still wrong. After editing and saving the edit in QGIS the shp appear fine in Arcgis and can now be normally edited.


texts: Is there a straightforward way to tweak LatLon on a sphere to represent true LatLon?

texts:
I have no experience in mapping, so please excuse me if this sounds foolish.

I'm trying to use Microsoft MapPoint (2004) software to geocode a number of addresses in Excel. However, my version of MapPoint does not have a nice, neat way to return the latitude and longitude of an address; it has to be calculated.

The method I've found to calculate latitude and longitude, in a nutshell, uses the great-circle distances from the location to the North Pole and the prime meridian. These distances are calculated using spherical coordinates, so they're fairly accurate, but they're not as precise as I'd like.

Is there a way to tweak these, let's call them "spherical" coordinates, so that they better represent true "elliptical" coordinates?


texts:
I have an exercice consisted on:

X is a binary with raw data from the Austrian radar services.
Each radar station scans 16 vertical surfaces with a displacement of 1 km.

The file contains the reflectivity data for all locations where radar stations are available.
This is called the radar composite for Austria.

I have already codified the binary file following this:

Byte           Example                Description
0              1                      Format id (1)
1-7            31-Mar-1998 11:55      Date/Time
18-19                                 Number of rows (short int LSB first) (=NY)
20-21                                 Number of columns (=NX)
22                                    Number of 3D-levels (=L)
23-24                                 distance(m) between the 3D-levels (short int LSB first)
25                                    Number of quantisation levels (=N)
26-26+n*2-1                           quantisation levels (N * short int LSB first) in 1/10mm/h
X1-Y1                                 Pixels (NX*NY Bytes) from upper-left to down-right on the 1st 3D-Level
X2-Y2                                 Pixels (NX*NY Bytes) from upper-left to down-right on the 2nd 3D-Level
XL-YL                                 Pixels (NX*NY Bytes) from upper-left to down-right on the Lth 3D-Level

The problems asks:

Create a picture (PNG) with the highest value for each column (vertically)

Do you understand what does it ask based on the data I have?
Value I assume comes from the Pixels value (which is codified in the range [1,14]) and map it into the quantisation_levels but I do not understand the question as a whole


texts: serving UTFgrid information made with mb-util with leaflet

texts:
I've created the mbtiles file with tilemill and exported it on my ubuntu webfolder with the following command:

 mb-util --scheme=xyz /folder/tile.mbtiles /var/www/html/my-tiles

To get utfgrid information in my leaflet map itself I used danzels utfgrid plugin and the following code which is a nice copy of danzels example:

var tiles = L.tileLayer('http://localhost/my-tiles/{z}/{x}/{y}.png', {
            maxZoom: 9,
            attribution: '&copy; OpenStreetMap contributors'
        });

        var utfGrid = new L.UtfGrid('http://localhost/my-tiles/{z}/{x}/{y}.grid.json');

        var map = L.map('map')
                .setView([50,10], 5)
                .addLayer(tiles)
                .addLayer(utfGrid);

A common error is written to the console:

ReferenceError: grid is not defined

So I've checke my grid.json files and they have the following coarse structure:

grid({"keys": ["", "6402"], "data": {"6402": {"namen": "cool name", "attrib 1": 0.123456}}, "grid": [" ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " !!! ", " !!!!! ", " !!!!!!! ", " !!!!!!!! ", " !!!!!!!!! ", " !!!! !!!! ", " !! !!!!! ", " !!! !!!!! ", " ! ! ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " "]});

Which looks fine for me.
So what am i missing? In the example from milkator he is using:

        var utfGrid = new L.UtfGrid('http://{s}.tiles.mapbox.com/v3/milkator.press_freedom/{z}/{x}/{y}.grid.json?callback={cb}')

Is there something to keep in mind using mb-util or anything else? My goal is to get the infromation from the UTF grid in a info box...


texts: Is there any way to use a regular expression to filter a string when using CartoCSS expressions?  I have geojson of some ways and route relations I obtained using Overpass Turbo.  When creating a new vector tile source from this geojson, Mapbox Studio flattens the nested @relations property to a string.  I'd like to use one of the relation's tags for labeling, but it appears I'll need a regular expression to extract it from the flattened string.  Is this possible, or are only simple expressions supported?  If not possible does anyone have suggestions on dealing with relations in CartoCSS (perhaps some imposm pre-processing is necessary but I'd like to avoid it if possible).


texts: Install QGIS Server and Web Client under Ubuntu 14.04

texts:
I tried to install QGIS Server and QGIS Web Client on my Ubuntu machine by following this tutorial http://www.aboutgis.ro/blog/wp-content/uploads/2013/04/qgis_desk_web_server_client_tutorial.pdf.
At the end of the tutorial i should start the browser with localhost/qgis-web-client/site/index.html. But i get the error that the /qgis-web... was not found.
the localhost works.
does anybody nows what's wrong woth the tutorial i followed step by step after: -->
install needed packages by entering in the terminal:
Edit1:

at the last step i copy the folder into /var/www/html/

now i can see the qgis-web-client/index.html but when i want to start the helloworld or the other map i get the following eror by firebug.
( TypeError: layerTree.root.firstChild is null )

it seems to me that the qgis server doesn't work correctly. I don't get any GetCapabilites request.

that's why i try to install qgis server from anitagraser http://anitagraser.com/2012/03/30/qgis-server-on-ubuntu-step-by-step/ but nothing changed because the files existing.


texts: Polygon overlay, calculating changes, in QGIS?

texts:
I'm trying to create a relationship file between two vintages of census tracts. (in this case, 2010 and 2013.)  Reason being to load ACS (American Community Survey) census data and use weighted means to tie it back to 2010 census tracts.  (The actual data for my project uses 2010 census tracts as of 2010 as a legal requirement.)

In truth, there's very little difference between these two census vintages (maybe 5-10 naming changes per year, and 1-2 geometry changes, in around 74,000 tracts.)  But I would rather not have to deal with these changes by hand, as more changes will be forthcoming as the ACS people continue to muddle with the geometries. I intend to use every year of ACS moving forward, and backward to 2010.

The shapefiles come from Census at the following source.  For the duration of this post, i'll be referring to FIPS, let me say that an 11 digit FIPS code uniquely identifies a census tract. (digits 1-2 identify state, 3-5 identify county within state, 6-11 identify census tract within county)

ftp://ftp.census.gov/geo/tiger/TIGER2010/TRACT/2010/

ftp://ftp.census.gov/geo/tiger/TIGER2013/TRACT/

Ideally, i want to be able to overlay them and get a sort of a table output, like:

FIPS_2010 | FIPS_2013 | AREA(Sq.KM)

Then quick SQL commands (area/sum(area) group by FIPS_year) would allow me to get:

FIPS_2010 | FIPS_2013 | AREA(Sq.KM) | PCT_2010 | PCT_2013

Then for any particular FIPS, if PCT_2010 = 1, then that element represents the complete tract for 2010.  If it's less than 1, then a geography change occured.

I'm just lost on how to do the initial overlay analysis in QGIS.


texts: EPSG 4326 to 3857 works on X axis but not Y with PROJ4NET

texts:
I am trying to convert a EPSG:4326 Lat Long into a EPSG:3857 reference. I have looked all around and have found slight differences in WKTs but nothing with different values.

The following code provides me with an accurate X coordinate according to here.

    Dim EPSG3857 As String = "PROJCS[""WGS 84 / Pseudo-Mercator"",GEOGCS[""WGS 84"",DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563,AUTHORITY[""EPSG"",""7030""]],AUTHORITY[""EPSG"",""6326""]],PRIMEM[""Greenwich"",0,AUTHORITY[""EPSG"",""8901""]],UNIT[""degree"",0.0174532925199433,AUTHORITY[""EPSG"",""9122""]],AUTHORITY[""EPSG"",""4326""]],PROJECTION[""Mercator_1SP""],PARAMETER[""Latitude_of_origin"", 0],PARAMETER[""central_meridian"",0],PARAMETER[""scale_factor"",1],PARAMETER[""false_easting"",0],PARAMETER[""false_northing"",0],UNIT[""metre"",1,AUTHORITY[""EPSG"",""9001""]],AXIS[""X"",EAST],AXIS[""Y"",NORTH],EXTENSION[""PROJ4"",""+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs""],AUTHORITY[""EPSG"",""3857""]]"

    Dim EPSG4326 As String = "GEOGCS[""WGS 84"",DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563,AUTHORITY[""EPSG"",""7030""]],AUTHORITY[""EPSG"",""6326""]],PRIMEM[""Greenwich"",0,AUTHORITY[""EPSG"",""8901""]],UNIT[""degree"",0.01745329251994328,AUTHORITY[""EPSG"",""9122""]],AUTHORITY[""EPSG"",""4326""]]"

    Dim srcCRS As ICoordinateSystem = Converters.WellKnownText.CoordinateSystemWktReader.Parse(EPSG4326)
    Dim tgtCRS As ICoordinateSystem = Converters.WellKnownText.CoordinateSystemWktReader.Parse(EPSG3857)

    Dim ctFac As New ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory
    Dim trans As Transformations.ICoordinateTransformation = ctFac.CreateFromCoordinateSystems(srcCRS, tgtCRS)

    Dim fromPT() As Double = {longitude, latitude}
    Dim toPT() As Double = trans.MathTransform.Transform(fromPT)

    Return toPT(0) & " " & toPT(1)

Using coordinates

longitude = -77.056610183150184, latitude = 38.871810818181821

it produces

-8577902.61, 4676494.62

instead of

-8577902.61, 4703326.17

Where am I going wrong?


texts: Why is Utility Network Analyst toolbar grayed out?

texts:
When I pull up my Utility Network Analyst toolbar I cannot use it because its grayed out. My method for pulling it up is Customize > Toolbars > Utility Network Analyst. I am not able to use any of the features of this toolbar.

Why is my Utility Network Analyst toolbar grayed out?


texts:
I'm currently working with a large number of CF-compliant NetCDF climate grids with global coverage that are organized as tiles across a file structure pertaining to particular "variables" and "years" (e.g. precipitation-year1950-001.nc, precipitation-year1961-125.nc, tmax-year2008-650.nc, etc.). They're currently organized in a sensible file structure (:project/variable/year/files.nc)

For each variable*year, I would like to spatially merge the tiles (similar in call to gdal_merge outputfile.tif /path/to/data/*.tif) using a walking method and I'm wondering what would be the most efficient/reliable way of doing so.

From my research, two (labour-intensive) methods come to mind:

Performing a walking loop in BASH that calls the Climate Data Operators (CDO) function mergegrid. Since mergegrid can only accept two input files for each output, even if I were to merge them efficiently, I would need to call it > 400x for each variable*year.
Performing a walking loop in Python with the NetCDF4 and NumPy libraries that creates an empty array the size of the the Earth, which reads each NetCDF to an np.array and writes it out to the new NetCDF, appends the metadata, and iterates through the variables*years.

Ideally, I need a method that is easily replicated after I've written the function out and if I can stay in Python3, that would be great. What I would like to know is which is the better way of merging lots of nested NetCDF data or is there another method/modules that I haven't come across that would be better at performing large spatial merges of NetCDF data.


texts: How to draw radius buffer around the points and perform intersection with other locations list in GeoPandas Python?


texts: Appending field from file geodatabase table to polygon feature class?

texts: I have an existing polygon feature class with about 30 attribute fields. I would like to append another 20 fields from an existing table in the same file geodatabase. Can this be done using an out-of-the-box tool? This is not a problem when creating a new feature class; in this case I would use the Import button on the Feature Class Properties:Fields tab.  Now that my feature class contains records (as well as a subtype field with subtypes, in case this makes a difference) the Import button is unavailable (dimmed).


texts: Using python to read each row and create one sided buffer depending on attribute table

texts:
I'm still fairly new to GIS and Python but I have a line feature class with one attribute called "View Side" which indicates which direction the buffer should be. I want to create a Python script to iterate through the attribute table and buffer according to the View Side. I have the following code so far, but the end results (the buffers) contain no attributes and I want all the attributes from the original line feature class to be preserved.

buffercursor = arcpy.SearchCursor(fc_albers)
for rows in buffercursor:
    if rows.getValue("View_Side") == "L":
        side = "LEFT"
    elif rows.getValue("Observation_side") == "R":
        side = "RIGHT"
arcpy.Buffer_analysis(rows.Shape, fc_buffers2, 250, side, "FLAT", "NONE")
arcpy.Append_management (fc_buffers2, target, 'NO_TEST', '', '')

Why am I doing this? This is a type of walk and view survey. I want to be able to append things to a "master file" down the line. The goal here is everytime I receive a shapefile, I could just plug it in and it'll create these buffers. I could always just select all that are L/R and buffer those and combine them, but this isn't what I'm looking for.


texts: Grab set of datapoints for each step for CartoDB torque timeSlider (through JavaScript)

texts:
So I have a regular CartoDB map that has a torque layer with a timeSlider. I know how to access the date for the slider object, but is there a way to access the datapoints for each step as it goes along?

What I want to try and do is access the set of datapoints for each 'step' of the timeSlider, which would be each day in this case. So as the slider moves or steps, I'd have access to the set of data points for that particular day/step.

window.onload = function () {
      cartodb.createVis('map', 'https://wkaravites.cartodb.com/api/v2/viz/e04e6a5e-8e2e-11e5-9b84-0e787de82d45/viz.json', {
          zoom: 11
        })
        .done(function(vis, layers) {
          var slider = vis.timeSlider;
        })
        .error(function(err) {
          console.log(err);
        });
    }


texts: How to specify two different and separate conditions in query builder?

texts:
For example, I have a table containing well log data for multiple operable units.  I would like to make a query outputting a table where Operable unit=001 OR 002, AND Formation= Residuum.  I would like to first select only values where operable unit equals 001 or 002, and then from those values select only those where the formation=residuum.

Currently the output from [OU] = '001' OR '002' AND [Formation] = 'Residuum', will output a table that will include other operable units as long as the formation=residuum.

I am using ArcMap 10.


texts: Why do tiles seems to jump around when using ArcGIS cache layer in Openlayers?

texts:
I am use Openlayers to add ArcGIS  tilecache layers using arcgiscache class or another way which is jsonp.

I can show the map properly, but when I pan, zoom or drag the map, the tiles seem to jump around, it does not look any good.

Can anybody give me some advice?
In this example it has been noticed that a problem the users encounter is that the tiles seem to "jump around", but I do not know how to handle it ！

var markers, layer;
var features = [];
var map, baseLayer,
    layerURL = "myurl";

function init() {
    var jsonp = new OpenLayers.Protocol.Script();
    jsonp.createRequest(layerURL, {
        f: 'json',
        pretty: 'true'
    }, initMap);

}
function initMap(layerInfo) {
    baseLayer = new OpenLayers.Layer.ArcGISCache("arciscache", layerURL, {
        layerInfo: layerInfo
    });
    map = new OpenLayers.Map('map', {
        maxExtent: baseLayer.maxExtent,
        units: baseLayer.units,
        resolutions: baseLayer.resolutions,
        numZoomLevels: baseLayer.numZoomLevels,
        tileSize: baseLayer.tileSize,
        eventListeners: { //set event handler
            "moveend": moveEvent}

        //displayProjection: baseLayer.displayProjection
    });

    var zoomControl = map.getControlsByClass('OpenLayers.Control.Zoom')[0];
    map.removeControl(zoomControl);

    map.addLayers([baseLayer]);
    map.addControl(new OpenLayers.Control.MousePosition());

    map.addControl(new OpenLayers.Control.Navigation());
    map.addControl(new OpenLayers.Control.PanZoomBar());
    if (!map.getCenter()) {
        map.zoomToMaxExtent();
    }
}


texts:
I'm a (reasonably) proficient GIS user (mostly QGIS, but some Arc GIS) who has been thrown the conundrum of converting/importing CAD urban planning drawings (2D drawings produce in AutoCAD 2014) to GIS shapefiles (or similar).

I have had some success in importing the drawings, but have realised that all of the drawings have been produced as lines (not polylines or polygons) which means I am losing a layer of detail that would be really very useful in GIS.

For example, I have a site plan of a new housing layout which shows hundreds of individual plots of land, which 'look' like polygons but are actually represented in AutoCAD as lines. As I understand it from reading through forum posts here it's probably best that these plots are converted to polygons in AutoCAD before we attempt to export/import them to GIS?

I'm really just looking for a solid enough workflow to produce 2D architectural plans in AutoCAD and import them to GIS so that I can now attach additional attributes to individual polygons (e.g. buildings, plots) and polylines (roads, rivers etc).

From first glance the CAD files are a bit of a mess (please see example DWG file here). They are also huge e.g. new housing estates with +4000 plots in some cases. I think refining the workflow at the CAD end is probably something that needs to happen. My other constraint is that with the exception of AutoCAD 2014 I'm pretty much restricted to using opensource software (QGIS etc)

Has anyone come across any useful resources that would point me in the right direction?!


texts: Oracle Maps filter.In allows me to only remove items from the map

texts:
Working from a basic example generated by Oracle MapBuilder. I have a base map and a FOI layer.
At the start no filter is applied and the map shows all object types (column OBJECTTYPE) available in the FOI layer.
 function addThemeBasedFOI()
  {
    layers["themebasedfoi1"] = layer1;
    map.addLayer(layer1);
    layer1.setZoomLevelRange(0,9);
   }

A checkbox, initially unchecked, calls this function:
  function setFilter(item)
  {
    if (item.checked) {
      layer1.applyFilter(new OM.filter.In('OBJECTTYPE',["WGD","PND","OTD","BTD"]),true);
    }
    else
    {
    layer1.applyFilter(new OM.filter.In('OBJECTTYPE',["WGD","PND","OTD"]),true);

    }
  }

First click on the checkbox checks it and all objecttypes are removed from the layer except the four specified ["WGD","PND","OTD","BTD"].
Second click on the checkbox unchecks it and the layer now show only the three object types. ["WGD","PND","OTD"]
With the checkbox checked again I expect the objecttype "BTD" to return but this doesn't happen.
How should this be done?


texts: Adding results from geoprocessing scripts to ArcGIS Pro?

texts:
I am in the process of updating scripts and models from ArcGIS Desktop to work in Pro however I'm having difficulty getting results from geoprocessing operations to display in the map window.

As an example, for the following snippet:

import arcpy
arcpy.env.overwiteOutput = True
aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.listMaps("*")[0]
inputFC = arcpy.GetParameterAsText(0)
outputWorkspace = arcpy.GetParameterAsText(1)
outputFC = arcpy.GetParameterAsText(2)
outputFC = outputWorkspace + "\\" + outputFC
arcpy.MultipartToSinglepart_management(inputFC, outputFC)
outputLyr = arcpy.MakeFeatureLayer_management(outputFC, "outputLyr")
m.addLayer(outputLyr, "TOP")

I get the error: "ValueError: outputLyr". The process completes successfully, it's just that I can't get the result to add to the map.

My knowledge of python is quite limited but as I understand it, a feature class needs to be made into a layer before being added to the map. Previously in Desktop we used the following:

layer_file = arcpy.mapping.Layer(outputFC)
arcpy.mapping.AddLayer(dataframe, layer_file)

Where am I going wrong?


texts: Simple Question: Raster from ArcGIS to Geoserver

texts:
we have a couple raster files which could be open in ArcGIS. They come with the extension of .aux and .rrd and they are the GRID file format from ESRI. We want to put them on Geoserver to overlay them on Google Map. They are all in WGS 84 projection.  I was wondering:

What is the flow chart of work we should do to finish this task?
What kind of conversions we need to do for the data?

We want relative fast loading of those data, so we need caching on the server. I would like to hear from you and try to write a tutorial thereafter.


texts: Problem creating viewshed raster in QGis - output black

texts:
I am trying to do a viewshed analysis in QGis using the viewshed analysis plugin in QGis repository (https://hub.qgis.org/projects/viewshed/).

I had as base elevation contours in shp format with elevation data in ELEV (1).

which I converted to a raster DEM (raster > analysis > grid), creating a raster layer (elevation_raster) (2).

Then I created a new point layer (for the observer) in the same coordinate system as the elevation raster (3).

Finally, I tried to generate a viewshed analysis with the given plugin, defining the elevation raster, the observation point, search radius (5000), observer height (7000), target height (9000), trying to output a binary viewshed as raster with comulative option (4).

The result (viewshed_comulative) is just black, and I dont know where to find logs / debug this situation (5).

Tried in both QGis 2.10 and 2.12.


texts:
I have created polygon graticule and I used ''polygons to lines'' to change the polygons into lines.
I want to change their azimuth to make them more feet into my polygons.
I had the idea to calculate the azimuth and then to find a tool to change the angle of the graticule according to the azimuth. I used this formula for the calculation-
degrees(azimuth(start_point($geometry), end_point($geometry)))

the problems are-

The calculated Azimuth is very weird, for most of the lines I got zero and in some I got numbers as you can see here-

seems like the lines didn't break very well, this is what happen when I use ''select'' and click on one line-

I don't find a tool that can help me to change the azimuth.

My end goal is to have those lines fit to my polygons as much as possible.if you have any other -more-simple idea of how to create this net inside my polygons. please share with me


texts: How to geocode addresses from MySQL and add them like markers to the Google map?

texts:
I'm using Google Maps API and Google Geocoding for geocoding addresses. I have table named 'users' and column named 'location' in MySQL database. Column 'location' contains addresses. I need to connect to the database, then geocode those addresses and add them like markers to the map.

First I connected to the database with PHP:

  <?php

    // connection to mysql

    $username="someusername";
    $password="*****";
    $database="GISwork";
    $hostname = "somehostname.com";
    $table = "users";

    $conn = mysql_connect($hostname,$username,$password) or die("Connecting to MySQL failed");

    mysql_select_db($database, $conn)
    or die("Selecting MySQL database failed");

    if (!mysql_select_db($database))
    die("Can't select database");

    // sending query
    $result = mysql_query("SELECT location FROM {$table}");
    if (!$result) {
        die("Query to show fields from table failed");
    }
    ?>

Then I defined geocoding and markers:

var userLocation = //it should be that $result from database

        var geocoder = new google.maps.Geocoder();
            geocoder.geocode( { 'address': userLocation}, function(results, status){
               if (status == google.maps.GeocoderStatus.OK) {
                  map.setCenter(results[0].geometry.location);
                  var marker = new google.maps.Marker({map: map, position: results[0].geometry.location });
               } else {
                  alert("Geocode was not successful for the following reason: " + status);
               }
            });

    }
    google.maps.event.addDomListener(window, 'load', initialize);

My question is:
a) How to pass php variable $result to the javascript variable userLocation to be readable for Google Geocoding?
b) Is that function for geocoding and adding markers enough for all markers or should it be different regards to I need to import all addresses from the database (multiple markers)?
c) How to add info windows to all markers?

Thank you.


texts: Understanding Google Earth Pro ruler tool algorithms?

texts:
When using Google Earth Pro ruler (Tools / Ruler or by clicking on the "ruler" icon), the ruler tool tells us the heading from one point to another which is not the same as that obtained by Karney's algorithm (a few tenths of degrees for long distances).

Could taking the geoid into account explain this difference?

Are there any algorithms to calculate the heading with the same precision as Google Earth Pro?


texts: Can I sell maps to third parties if I have Google Earth Pro License?

texts:
Not sure if this question is suitable for SE but I can't find answers anywhere else.

I work as a GIS Analyst for a mining company that has its own GEPro license.  As part of my job, I make PDF maps, which are then used externally on the company's website for promotional purposes.  The maps display geological information, with GE satellite imagery as a base layer.

Per my understanding of GEPro's TOS, this is acceptable.  I am an employee of a company which has a valid license, and I am making maps on their behalf.

Now for the tricky part.  I am starting to do freelance/contract map work for other companies.  I have my own GEPro license, and would be delivering to my clients PDF maps of geological data with satellite imagery as the base layer.  The clients would then use the PDF maps externally, for promotional purposes. I think that this would classify my clients as "Third Parties".

As far as I see it there are three options:

1) Making maps for third parties for a fee is acceptable with a GEPro license

2) Making maps for third parties for a fee is acceptable with a GE Enterprise license

3) I can make as many maps as I want, so long as the company I'm making them for has a valid GEPro license

...

I've read over most everything Google has put out as far as TOS goes, but clarity is hard to come by.  The GE Enterprise TOS talks about not selling "Intermediate Data Products" or "Fused Databases", but I dont think those are applicable to what I would be providing (though admittedly, I am not crystal clear on what those two terms mean).

Anyone have any insight/experience with this?

Thanks for reading


texts: qgis2web export categorized renderer symbology

texts:
I'm trying to export my QGIS (2.14.5) categorized renderer symbology to a leaflet project with qgis2web (1.19.0).

My example QGIS project looks like this:

The resulting (preview, respectivly) leaflet looks like shown below. Note that the legend shows categories, but not the map:

qgis2web help mentions only that in leaflet only the first symbol layer is displayed, but that is not the issue in this case, me thinks.

Do I miss something, or do I something wrong?


texts: QgsProcessingException: There were errors executing the algorithm

texts:
I have been trying to run the nearest distance to hub points process (https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/vectoranalysis.html#qgisdistancetonearesthubpoints) in Jupyter notebook.
As of the data, I have 2 geojson files (polygon) i.e. inputvector2 and inputvector3 in this case. I have attached the error I get as image below.                        PS: I am a newbie

Code snippet
vertice_out = r"D:\nearesthub.geojson"
# just delete the output if exist
if os.path.exists(vertice_out):
    os.remove(vertice_out)
else:
    pass
params = {
    'INPUT': inputvector2,
    'HUBS' : inputvector3,
    'FIELD' : 'name',
    'UNIT' : 0,
    'OUTPUT': vertice_out
}
feedback = QgsProcessingFeedback()
processing.run("qgis:distancetonearesthubpoints", params, feedback=feedback)


texts: Buffering in pyQGIS? Task: Find all countries that are in the range of missiles

texts:
I am beginning user of pyQGIS and I have task to do. I need to find all cities outside country called R, which are in range of missiles based in city K. The range of missiles is 250 km. Missiles located 15 km inside the boundary of the city K.
I loaded Layers:

I. Cities:

layer = QgsVectorLayer("C:\cities.shp","miasta", "ogr")
layer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(layer)

II. Countries:

layer2 = QgsVectorLayer("C:\country.shp","panstwa","ogr")
layer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(layer2)

Then, I reaveled all cities outside R:

for f in layer.getFeatures():
geom = f.geometry()
if f['CNTRYNAME'] != 'R':
    print f['NAME']

And their coordinates:

for f in layer.getFeatures():
geom = f.geometry()
if f['CNTRYNAME'] != 'R':
    print geom.asPoint()

Now I don't know how to find out how many cities (except ones in R) are in the range of missiles.

How to create buffer which is 15 km from the border and how can i find all cities 250 km from that buffer? Maybe there is other way of doing that task?

I also have to find time that the missile travels from boarder of Kaliningrad to Lublin, Poland and ellipsoidal distance. Missile's speed = 2.5 Ma (850 m/s)
I have no idea what to do next.


texts: Create Statewide Hydrologic Soil Group Map

texts:
This may be better for Engineering? Not sure.

Background:
Previously I have downloaded SSURGO soil data for the state of Pennsylvania by downloading each county and merging the shapefiles. I then combined the mapunit and muaggatt tables and joined them with the spatial data based on the soil name (e.g. "GbC", "WhB"). I use this to create thematic maps of Farmland classification and Hydric conditions programmatically using arcpy.

Goal:
Now I would like to do the same to create thematic maps of Hydrologic Soil Classification. It looks like the appropriate table is comp.txt however this table doesn't appear to have a similar soil name field.

Question:
Has anyone done this before? That is, linked SSURGO soil spatial data to appropriate tables to display Hydrologic Soil Classification? Or does anyone have guidance on any pathways that would be best to look into?

As pointed out in the answer to the question linked below, some have used PostgreSQL to connect to GIS, (would be fantastic) but I have no experience using databases and I can't find any descriptions of how exactly to do this. Can I download the entire NRCS SSURGO database?


texts: How to georereference a TIFF file which has already a spatial reference?

texts:
I have a TIFF file of Poland and a self drawn vector map of Poland. The vector map has no spatial information.

I am intending to reference the TIFF exactly to the shape of the self drawn map.

My problem is, that the TIFF is not choosable in the georeferencing tool. I suspect it is because it already has a CRS and so on. Unreferenced data like a JPEG can be choosed and referenced.

Is there a way to rereference my raster data?

It is very important to contain the exact shape of the self drawn map and I think the geoeferencing tool should do well to fix my raster appropriately. It is furthermore important to contain the classification info, it would not be suficient to convert it to a jpg.


texts: Mapinfo SQL - create table that shows where a new polygon does not intersect with another

texts:
I am attempting to carry out SQL to - create a table that shows where a new polygon (Proposed BMO Overlays) does not intersect with another (Existing BMO Overlays).

I have attached an image on my Map layout to help show the layers involved.
Can you please help?


texts: Finding and replacing with if then else statement in field calculator of ArcGIS Pro?

texts:
I have a column of strings that always start with the first two characters being numbers, like 39 or 07, etc.  In a new column I want to process the original column so that I replace any of those first two numbers that are not 39 with a 39, but making sure the remaining trailing string characters are still there. So a find and replace but index on the first two characters.  I want to use the .replace() command but not sure how to tie in the index, replace command and an If Than Else statements to make it work.

In field calculator it would be like this:

replace(!Permit_num!)

code block:

import re
def replace(val):
    if (val[0:1] == "07"):
        return re.sub ('07','39'+val[2:], val)
    elif (val[0:1] == "10"):
        return re.sub ('10','39'+val[2:], val)
    elif (val[0:1] == "11"):
        return re.sub ('11','39'+val[2:], val)
    elif (val[0:1] == "19"):
        return re.sub ('19','39'+val[2:], val)
    elif (val[0:1] == "24"):
        return re.sub ('24','39'+val[2:], val)
    elif (val[0:1] == "27"):
        return re.sub ('27','39'+val[2:], val)
    elif (val[0:1] == "34"):
        return re.sub ('34','39'+val[2:], val)
    elif (val[0:1] == "35"):
        return re.sub ('35','39'+val[2:], val)
    elif (val[0:1] == "39"):
        return re.sub ('39','39'+val[2:], val)
    elif (val[0:1] == "50"):
        return re.sub ('50','39'+val[2:], val)
    elif (val[0:1] == "54"):
        return re.sub ('54','39'+val[2:], val)
    elif (val[0:1] == "56"):
        return re.sub ('56','39'+val[2:], val)
    elif (val[0:1] == "99"):
        return re.sub ('99','39'+val[2:], val)

I tried the above, but does not work, get a 99999 error and says my field I am calculating into is not nullable.  The new field is also a string.


texts: shp2pgsql no more AddGeometryColumn() function?

texts:
Problem:

I encounter a fresh error when using shp2pgsql on Ubuntu 18.04.
This commands worked a few months ago:

$ sudo shp2pgsql -s 4326 /path/to/myshapefile.shp | psql -h localhost -p 5432 -d databasename -U postgres
Field height is an FTDouble with width 31 and precision 15
Shapefile type: PointZ
Postgis type: POINT[4]
Password for user postgres: ************************************
SET
SET
BEGIN
CREATE TABLE
ALTER TABLE
ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 1: SELECT AddGeometryColumn('','myshapefile','geom','4326'...
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
(...)

Question:

Any idea on how to fix this?

Environment:

Ubuntu 18.04 64 bits: 5.0.0-37-generic x86_64 GNU/Linux

SELECT Version();
"PostgreSQL 10.10 (Ubuntu 10.10-0ubuntu0.18.04.1) on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0, 64-bit"

SELECT postgis_version();
"2.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"


texts: Any consideration of moving QGIS project file to SQLITE form XML?

texts: Just curious, has anyone ever proposed storing the ".qgs" project file in SQLITE instead of XML?  Seems like it might be faster, easier to parse, and you could hide the changes behind an API.  But I know I am not volunteering... So, if anyone has heard such a discussion, please post the links or comment.


texts: Convert a NetCDF layer to GeoTIFF image using gdal_translate, applying mask from another layer

texts:
My objective is to convert a layer in a NetCDF file to GeoTIFF image with gdal_translate, applying mask from another layer in the same NetCDF file.
Sample NetCDF file: sample file on dropbox (~ 6Mb)
GDAL version:  3.1.1, released 2020/06/22. Windows OS: 10
Layers in the test.nc file using gdalinfo:
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"test.nc":quality_level
  SUBDATASET_1_DESC=[1x9000x18000] quality_level (8-bit integer)
  SUBDATASET_2_NAME=NETCDF:"test.nc":l2p_flags
  SUBDATASET_2_DESC=[1x9000x18000] l2p_flags (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"test.nc":or_number_of_pixels
  SUBDATASET_3_DESC=[1x9000x18000] or_number_of_pixels (16-bit integer)
  SUBDATASET_4_NAME=NETCDF:"test.nc":sea_surface_temperature
  SUBDATASET_4_DESC=[1x9000x18000] sea_surface_subskin_temperature (16-bit integer)
  SUBDATASET_5_NAME=NETCDF:"test.nc":dt_analysis
  SUBDATASET_5_DESC=[1x9000x18000] dt_analysis (8-bit integer)
  SUBDATASET_6_NAME=NETCDF:"test.nc":satellite_zenith_angle
  SUBDATASET_6_DESC=[1x9000x18000] satellite_zenith_angle (8-bit integer)
  SUBDATASET_7_NAME=NETCDF:"test.nc":sses_bias
  SUBDATASET_7_DESC=[1x9000x18000] sses_bias (8-bit integer)
  SUBDATASET_8_NAME=NETCDF:"test.nc":sses_standard_deviation
  SUBDATASET_8_DESC=[1x9000x18000] sses_standard_deviation (8-bit integer)
  SUBDATASET_9_NAME=NETCDF:"test.nc":wind_speed
  SUBDATASET_9_DESC=[1x9000x18000] wind_speed (8-bit integer)
  SUBDATASET_10_NAME=NETCDF:"test.nc":sst_dtime
  SUBDATASET_10_DESC=[1x9000x18000] sst_dtime (32-bit integer)

Convert the band with name "sea_surface_temperature" to geoTiff:
gdal_translate -of GTiff NETCDF:test.nc:sea_surface_temperature test.gtiff

The operation is straightforward up to this point. But, my requirement is to apply the bit-mask in the layer "l2p_flags", specifically the 5th bit which specifies whether a pixel is in Day or in Night.
How can one specify in 'gdal_translate' to use the 'l2p_flags' layer as mask and apply the 5th bit?


texts: gdal_translate changes colors converting qgis generated PDF to GTiff

texts:
I am producing a geospatial PDF using Qgis 3.10.4. One of my polygon layers is symbolized using red diagonal hatching. The export to pdf gives the correct rendering -

I then uploaded the pdf to the Avenza Map store, and the colors of the hatching changed to cyan or something close. Wondering whether there might be some CMYK vs RGB thing going on (though that makes no sense since no other colors changed), I tried converting the pdf to GTiff using gdal_translate:

gdal_translate -of GTiff --config GDAL_PDF_DPI 300 "Llaves NM.pdf" Laves_color_Test.tif

The result looks like this:

So it's not just Avenza's rendering stream that's changing the colors.

To complete the tests, I tried exporting directly to .tif and get this, which is correct:

I would much prefer to use pdfs, but not if the colors are going to shift.
Is there a setting or something that fixes this issue?


texts: How to cluster points of interest by lat lon for my WP7 app?

texts:
I have alread asked this question on SO, so I apologize if duplicating this is bad form. I will link to each of the sites if I get the right answer.

I have a collection of my own PointOfInterest classes (approx. 1500) each having their own Latitude and Longitude double properties.

Im trying to draw them on my screen, but at a certain logical zoom level, there is no point showing some of them because they are so close together.

How can I very efficiently group POI's by their lat lon properties?

I have this type of thing:

var pointOfInterests = (from p in PointsOfInterest select p).Distinct(new EqualityComparer()).ToList();

where the EqualityComparer is:

public class EqualityComparer : IEqualityComparer<PointOfInterest>
{
    public bool Equals(PointOfInterest x, PointOfInterest y)
    {
        return Math.Round(x.Latitude.Value, PointOfInterest.DecimalPlaceFilterLevel) == Math.Round(y.Latitude.Value, PointOfInterest.DecimalPlaceFilterLevel) &&
            Math.Round(x.Longitude.Value, PointOfInterest.DecimalPlaceFilterLevel) == Math.Round(y.Longitude.Value, PointOfInterest.DecimalPlaceFilterLevel);
    }

    public int GetHashCode(PointOfInterest obj)
    {
        return Math.Round(obj.Latitude.Value, PointOfInterest.DecimalPlaceFilterLevel).GetHashCode() ^ Math.Round(obj.Longitude.Value, PointOfInterest.DecimalPlaceFilterLevel).GetHashCode();
    }
}

and the PointOfInterest.DecimalPlaceFilterLevel is a static int property that I set when the user is at a certian zoom level.

But this isnt working, I keep getting overlapping POI's and its not very fast... since I am on the phone, I need it to perform very well.


texts:
I want to connect a PostGIS database to ArcMap 10.4. I work on Windows 7 with a PostgreSQL 9.6. Interoperability extension is installed and activated. I can see the PostGIS layers in the Catalog. I can see the field names of the layers. Even the preview works:

And also in QGIS it works fine.
But whenever I try to load one of these PostGIS layer into my project I get an error message: "Attribute column not found" and ArcMap does not show the feature geometries and the attribute columns:


texts: PostGIS Split Multipolygon with another Multipolygon

texts:
Since PostGIS does not have the ability to split multipolygons by other multipolygons or multilinestrings I found this workaround

CREATE TEMP TABLE boundaries AS (
SELECT
  a.statefp as state,
  a.countyfp as countyfp,
  a.name as county,
  a.cd110fp as cd,
  CASE WHEN ST_Within(a.geom, b.geom) THEN a.geom
  ELSE ST_Multi(ST_Intersection(a.geom, b.geom)) END AS geom
 FROM us_cd2008 b
JOIN countyparts2008_with_cd a ON ST_Intersects(a.geom, b.geom));

CREATE TEMP TABLE boundaries2 AS (
SELECT * FROM (
SELECT state, cd, countyfp, county, (ST_DUMP(b.geom)).geom FROM boundaries b) as temp
   WHERE ST_GeometryType(temp.geom) = 'ST_Polygon');

CREATE TABLE countyparts2008_split_v2 AS (
SELECT ST_Multi(ST_Union(geom)) AS geom,
     state,
     cd,
     countyfp,
     county
   FROM boundaries2
   GROUP BY cd, countyfp, county, state);

However the ouput looks like this...

Not sure how to fix it so that postgis doesn't include those lines towards the bottom. Any idea?


texts: How to correctly use Datepart( ) in ArcGIS?

texts:
I have two date fields dd.mm.yyyy with a start and an end date in my attribute table. Now I want only to extract the years from these fields to new fields. So I created another field, and used year_start_date = DatePart ( "YYYY", [Start_date] ) in the field calculator.

But the result is a complete date (dd.mm.yyyy) more than 100 years ago.

Just to test, I created a new field and added with the function Now() the actual date (13.05.2014) and tried again to use year_start_date = DatePart ( "YYYY", [Start_date] ) and I got "06.07.1905" as a result.

What am I doing wrong? Is there another way to do this?
I use ArcGIS 10.2.


texts: Counting geometry intersections between two linestrings

texts:
I am running PostgreSQL 11.2 and PostGIS 2.5.3
I have two linstrings intersecting each other like in the picture below:

The red line passes over itself and have vertices on, in order, A-B-C-B-D, so technically the red line intersects the black one twice, once from B to C and an second time from C to B.
I need to count how many times the red line intersects the black one.
I try ST_Numgeometries(ST_Intersection(red.geom, black.geom)) but only get 1 as a result. And obviously the ST_Intersection(red.geom, black.geom) returns only one geometry.
I can imagine convoluted ways to get the result I want, like for example decomposing the red line in simple lines AB, BC, CB and BD, and then run ST_Intersection, but since I want to run that in a spatial view I believe it will greatly slow down the process. Is there a simple way to achieve what I need, like a PostGIS function I didn't notice?
EDIT 1
Following the answer from Taras, here is an SQL example using his second option, which still returns 1 as a result :
WITH line1(geom) AS (VALUES ('LINESTRING(0 0, 0 1, 3 1, 0 1, 0 3)')),
line2(geom) AS (VALUES ('LINESTRING(1 0, 1 5)'))
SELECT SUM(ST_NumGeometries(ST_Intersection(l1.geom, l2.geom)))
FROM line1 AS l1
JOIN line2 AS l2 ON ST_Intersects(l1.geom, l2.geom)
GROUP BY st_intersection(l1.geom, l2.geom)


texts: Why QGIS cannot georeference a Grib2 file with CRS (PROJ.4 and OGC WKT2)?

texts:
Situation:

I have some .grib2 files that I downloaded from a (supposed to be) reliable source. When I do

$ gdalsrsinfo file.grib2

I get:

PROJ.4 : +proj=lcc +lat_0=38.5 +lon_0=262.5 +lat_1=38.5 +lat_2=38.5 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["unnamed",
  BASEGEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6371229,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]]],
  CONVERSION["unnamed",
    METHOD["Lambert Conic Conformal (2SP)",
        ID["EPSG",9802]],
    PARAMETER["Latitude of false origin",38.5,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8821]],
    PARAMETER["Longitude of false origin",262.5,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8822]],
    PARAMETER["Latitude of 1st standard parallel",38.5,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8823]],
    PARAMETER["Latitude of 2nd standard parallel",38.5,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8824]],
    PARAMETER["Easting at false origin",0,
        LENGTHUNIT["Metre",1],
        ID["EPSG",8826]],
    PARAMETER["Northing at false origin",0,
        LENGTHUNIT["Metre",1],
        ID["EPSG",8827]]],
  CS[Cartesian,2],
    AXIS["easting",east,
        ORDER[1],
        LENGTHUNIT["Metre",1]],
    AXIS["northing",north,
        ORDER[2],
        LENGTHUNIT["Metre",1]]]

So, from my understanding, the file is provided in a custom CRS, which is defined in two formats (PROJ.4 and in OGC WKT2) as part of the metadata of this file. (is that correct?)

Expectation:

According to the CF MetaData page, OGC WKT2 is used by GDAL and PROJ.4 is part of the OSGeo family.

So, as I have my QGIS 3.10 installed with OSGeo4W (Windows 10), I was expecting that, when I load this grib2 file into a QGIS, the layer would be loaded with a custom projection (but well projected) following the OGC WKT2 / PROJ.4 specifications.

Result:

QGIS says that this file has no CRS set and the layer is presented totally distorted.

Questions:

1-) Can I reproject this file to some other CRS (as ESPG:4326)?

2-) Do I need to install Proj4 or OGC WKT libraries separately on my OSGeo4W ?


texts:
I have a table that looks like the following

ID    Address                      Date
1     11 abc road                  19/03/2018
2     11 def road                  20/03/2018
2     11 def road                  21/03/2018
3     11 qaz road                  25/03/2018
4     11 zxc road                  26/03/2018
4     11 zxc road                  27/03/2018
6     11 qwe road                  18/03/2018
6     11 qwe road                  19/03/2018
6     11 qwe road                  20/03/2018

What I would like to do, is to merge all the rows that have the same id, into one record but merge the differing dates into that same record aswell.

So the final table would look like

ID    Address                Date
1     11 abc road            19/03/2018
2     11 def road            20/03/2018, 21/03/2018
3     11 qaz road            19/03/2018
4     11 zxc road            26/03/2018, 27/03/2018
6     11 qwe road            18/03/2018, 19/03/2018, 20/03/2018

What is the best pythonic way of achieving this?


texts: Adding table to ArcGIS MapServer so layers and tables can be associated?

texts:
I'd like to add tables to my MapServer. I want it to include both feature layers and tables which enables me to associate the layers with the relevant tables and retrieving the necessary fields of the tables. The resulting structure will be like the one in the following link;
http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Petroleum/KSPetro/MapServer

After clicking on the "tops" table, the relationship between the layers and the table "tops" can be seen. Though all the necessary tables and layers can be seen through ArcCatalog, it cannot be seen in the MapServer. After updating the ArcCatalog, The MapServer is reset.

Can anyone help me?


texts: OSM zoom levels correspond to what level of text labels?

texts:
I'm using Maperitive to render OSM map data of Mexico City as SVG files, destination Illustrator and then print (on paper, in a book).

At different zoom levels, the resulting maps have different levels of text.
At zoom=16 it looks like I get everything down to the last little street name, at 13, far fewer labels.

Since we want to print city maps with just the major categories - like names of suburbs - does anyone know what labels are included with the different zoom levels? If I knew that, I could do the maps without the thousands of street names that would be too tiny to read in print. OSM has something on it:

http://wiki.openstreetmap.org/wiki/Zoom_levels

...but not much.

Once I've worked it out, I want to do other maps of other cities, hoping that the same rules apply in all of them.


texts: QGIS and topological complexities/erros in polygons

texts:
I am using QGIS 2.14.19-Essen on Windows 10. I have a polygon shapefile of ecorregions in the Southern Cone of South America (Fig. 1), which is composed of 317 different features. I need to make it look like one whole polygon without internal limits. Then, I dissolved it using SAGA's Polygon Dissolve tool, which succeeded to some extet. As can be seen in Fig. 2, two lines remain. Line 1, which corresponds to the northern limit of Argentina, and line 2 which represents the borders between Argentina/Uruguay and Rio Grande do Sul in Brasil.

Line 1 can be easily eliminated using the node tool and selecting all the vertices to be deleted until none remains. Lines surrounding Rio Grande do Sul (No. 2) are a different story (Figures 3 and blown up detail in Figure 4) because the intersecting, zigzagging lines delimiting Southern Brasil (as a separate part) seem to preclude a successful dissolution.

Is there a straightforward way to get rid of those lines and nodes and have one single polygon? I've tried several alternatives using a wide range of pyhton/grass plugins, all to no avail.


texts: How do I calculate the distance among all points in set in QGIS?

texts: I am a very new user of any GIS.  In QGIS I have input my UTM data as a CSV file.  Under Vector -> Analysis Tools -> Distance Matrix, I have been able to create a distance matrix but it just takes one point and calculates the distance to each other point.  I want the distance between all possible points. Is there a simple way to do this?


texts: What typeface is used for County Boroughs on OS 6" to Mile Maps?

texts:
I came across this wonderful typeface on an OS 6" to Mile map (from National Library of Scotland):

By perusing the Characteristics Sheet it looks like it is representing a county borough (or country burgh in Scotland).

Does anyone know what the typeface is? and is it available to download anywhere? Or do I need to work out the minimal number of boroughs to get a full alphabet?


texts: GEE cloud-free Sentinel2 and linear Regression

texts:
GEE, Sentinel and linear Regression

I am using GEE (Google Earth Engine) to compute a Sentinel-2 based change analysis. What I would like is basically a linear regression per pixel, with the slope as an output band. However, I need to address issues with atmospheric contamination of imagery. Hence, I am running an image classification to then produce a mask based on my classes for atmospheric contamination (not happy with QA60).
I then tried to use

a) the linearFit reducer

Problem: If data is masked in only one scene within the collection, the linearFit is masked too (even though there might be other imagery available without mask)! However, I would like the masked data simply to be ignored as Null and the trend to be displayed using the rest of the available data. Hence I tried b), where beta can be specified as NULL

b) the rubustlinearRegression

Problem: I do not understand the output Layer! what is the weird crosses? If I use the pixel inspector it shows the blank regions as masked. Why is it masked differently? Is it because of the min and max values I specified for the lrImage layer? If that is the case, Is there a way to see the range of values within a band as to then define min and max correctly?
Generally, I am assuming that I cannot map the trend with clouds as noData because the cloud mask is run through the entire collection, which includes the time band. Would it be a solution to exclude the time band from the mask? If yes, how would I do that?

This is the link: https://code.earthengine.google.com/1cbaf9ad6265c166bdbe508c62d7ad4b

This is the code. It includes a point geometry = roi5, Sentinel2 collection = MSI, NSW = Polygon used as boundary, and Feature Collections as landcover classes.

//Compute a trend for the cloudfreecollection

// Select the bands to model with the independent variable first.
var trend = QAclclear.select(['system:time_start', 'NDVI'])
  // Compute the linear trend over time.
  .reduce(ee.Reducer.linearFit());

//the max: [-1, 1, 10] impacts the masking!!!!
Map.addLayer(
    trend,
    {min:0 , max: [-1, 1, 10000] , bands: ['scale', 'scale', 'offset']},
    'NDVI trend');

//Robust Linear Regression where beta (here Null) is the error outlier, and a default

// This function adds a constant band to the image.
var createConstantBand = function(image) {
  return ee.Image(1).addBands(image);
};

//Map the new constant band over the QAclclear (atmosph. dist. masked) collection, and select the constant, independent, and dependent value
var ConstantQAclclear = QAclclear
              .map(createConstantBand)
              // Select the predictors and the responses.
  .select(['constant', 'system:time_start', 'NDVI']);

print(ConstantQAclclear);

//Take the Collection with the constant value and apply the reducer
var RLR = ConstantQAclclear.select(['system:time_start', 'NDVI', 'constant'])
  // Compute the linear trend over time.
  .reduce(ee.Reducer.robustLinearRegression({numX: 2, numY: 1}));

print(RLR);

// The results are array images that must be flattened for display.
// These lists label the information along each axis of the arrays.
var bandNames = [['constant', 'time'], // 0-axis variation.
                 ['NDVI']]; // 1-axis variation.

// Flatten the array images to get multi-band images according to the labels.
var lrImage = RLR.select(['coefficients']).arrayFlatten(bandNames);

// Display the OLS results.
Map.addLayer(lrImage,
  {min: 0, max: [1, 1], bands: ['constant_NDVI', 'time_NDVI']}, 'OLS');

I am new to this.


texts: WMS request causes Geoserver Lat/Lon Bounding Box to Change

texts:
I have a layer with the following bounding box:

A bot (baidu spider) has been making the following request periodically:
/geoserver/si-x/wms/kml?layers=si-x:zabelli_leaf_ncep

This request causes the bounding box to change to the following:

Any ideas why? If I click "Compute from native bounds" and save the layer it will be fine until the bot hits the layer again at which point the bounds revert to the numbers in the second screen shot. If I make a similar request from a browser

http://myGeoserver/geoserver/wms/kml?&layers=gdd:agdd_50f

I don't receive an error, a .kml file downloads, and the bounds of the layer are messed up again. Any ideas on why this request would cause the Lat/Lon Bounding Box of the layer to change?


texts: Leaflet js dealing with a lot of markers in a single co-ordinate

texts:
I am using leafletjs and leafletjs market clustering.

I have a list of friends and I want to display them all on a map. The problem is that some of them life in a same house, so the co-ordinate for multiple markers would be the same. Issue occurs when there is more than 50 friends living in the same place.

Is there any way, that the markers could be hidden and when a cluster is clicked it instead would display a table containing all of the marker names?

My code that is adding markers:

export function markersFromData(map, markers) {
  return (data) => {
    const markerList = [];
    data.map((v) => {
      const title = v.name;
      const marker = L.marker(new L.LatLng(v.latitude, v.longitude), {
        opacity: 0,
      });
      marker.bindPopup(title);
      markerList.push(marker);
      return markers.addLayer(marker);
    });
    map.addLayer(markers);
    // eslint-disable-next-line
    const group = new L.featureGroup(markerList);
    map.fitBounds(group.getBounds());
  };
}

Map example


texts:
I´ve got the "OS Terrain 50" contours/spot-heights dataset.
I would like to output a series of maps with everything above 0m, 50m, 100m etc. from a region.
I know almost nothing about GIS, but have a reasonable technical knowledge.

There are lots of GIS programs, but which would be best/easiest to do this in.

Thanks

Edit: Limited (tagged) the program to the recommended QGIS.
I am wanting to make a physical model of some land. Terrain outline at different heights above sea level will give some templates to cut out and stack to make an actual model. Hope it´s not going too far off topic.


texts: Compute area coverage density estimate from very large polygon dataset

texts:
I am trying to come up with a pipeline for turning very large OSM PBF files (eg buildings layer for a whole country: ~100M polygons) into a more manageable density estimate (with arbitrary precision).
So far my approach has been to convert the PBF file to a more usable GPKG format with ogr2ogr:
ogr2ogr -select building -overwrite -progress -gt 5000000 -f GPKG buildings.gpkg buildings.osm.pbf multipolygons

followed by a geopandas script that loads chunks (in parallel) and sums up polygon areas for each tile (using polygon boundary centre as a quick tile assignment estimate):
precision_scale = 10

def parse_slice(s):
    gdf = gpd.read_file(gpkg_path, rows=s, driver='GPKG', mode='r')
    snapped = gdf.bounds.apply(lambda b: (
        int((b.minx+b.maxx)/2 * precision_scale),
        int((b.miny+b.maxy)/2 * precision_scale)
    ), axis=1)
    df = pd.DataFrame({'area': gdf.to_crs(CRS.from_epsg('6933')).area,
                       'snapped': snapped}).groupby('snapped').sum().reset_index()
    return df

merged_df = pd.DataFrame()
with multiprocess.Pool(processes = 4) as pool:
    for df in pool.imap(parse_slice, slices):
        merged_df = merged_df.append(df)

df = merged_df.groupby('snapped').sum().reset_index()

grid = gpd.GeoDataFrame(df.area, geometry=df.snapped.apply(
    lambda pt: Point(pt[0]/precision_scale, pt[1]/precision_scale)
), crs=gdf_crs)

This works fairly well for medium-sized data (PBF file < 2GB), but ogr2ogr seems to choke on anything bigger (writes up to ~9GB of GPKG file, then hangs forever without any error message).
Is there any method/tool/algorithm that is more cpu and/or memory efficient that I could use to achieve a similar result?


texts: Polyline selection of features with angle turns up to a specified angle

texts: I'm looking for an automated way to select line features that have specified angle between vertexes.  That in order to check water mains layer -steel pipes can have only 45 or 90 degrees turns and I need to check the digitizer work. ArcMap Data Reviewer have a tool for that, but i don't have this extension.


texts: Exporting mosaic dataset to geodatabase with arcpy

texts:
How would one go about exporting a Mosaic dataset from a enterprise geodatabase to a file geodatabase using ArcPy?

I'm able to copy a mosaic dataset via the ArcCatalog GUI by doing a copy/paste from one geodatabase to another, but I was hoping I could automate this with ArcPy.

I tried using the copy tool to do this but I can't seem to copy into an existing geodatabase:

import arcpy
arcpy.env.workspace = "<path to sde file>"
dest_gdb = "<path to new file gdb>" + "test.gdb"
arcpy.CreateFileGDB_management("<path to file new file gdb>", "test.gdb")
mds_list = arcpy.ListDatasets("*", "Mosaic")
    for ds in mds_list:
        arcpy.Copy_management(ds, dest_gdb)

The above snippet errors out with:

Traceback (most recent call last):
  File "migrate.py", line 85, in (module)
    arcpy.Copy_management(ds, dest_gdb)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py" line 3838, in Copy
  raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000725: Output data element: Dataset <path to file gdb>\test.gdb already exists.
Failed to execute (Copy).

I was also thinking of using the create mosaic dataset function and then  inserting the individual rasters from the source dataset into the new mosaic dataset, but I couldn't find any functions to list the rasters of a mosaic dataset.


texts: Creating polygons or map book from DBF file through ArcPy?

texts:
I have a dbf file, which contained a list of coordinate information as shown in below. It also could be thought as a catalog file. Example:

ID   XMIN       YMIN      XMAX       YMAX
1  -52800.0   -63000.0  -52000.0   -62400.0
2  -48800.0   -60600.0  -48000.0   -60000.0
3  -48800.0   -61200.0  -48000.0   -60600.0
4  -49600.0   -61800.0  -48800.0   -61200.0
5  -48800.0   -61800.0  -48000.0   -61200.0

So, how can i extract these information to create shapes or a map book as the picture shown in below. Using Python.


texts: QGIS and ogr2ogr not using distortion grid for NZMG by default

texts: is there a way to get QGIS and its projection engine OGR2OGR (GDAL software) to use a distortion grid by default. We are continually using a combination of NZMG and NZTM (NZGD2000 datum) datasets that we want to display and edit. We have noticed QGIS (ver 2.14) shows the same point on both projections to be misaligned by a few metres.  It appears ogr2ogr is using the 7 parameter algorithm as the default. We want it to use the nzgd2kgrid0005.gsb distortion grid (which it has in its ..\share\proj folder) by default. I didn't want to have to set up a Custom CRS if I could help it.


texts: How can I get a remote village transportation network from openstreetmap?

texts:
How can I get such a remote village transportation network from openstreetmap?

nodes represent remote villages
links represents public transportation tools, like buses

Something like that:

It is reasonable to regard bus stations as villages or towns. But how can I extract those info?

Should I create a map firstly as described in Openptmap/Installation (it is too complicate for me) and then filter public transportations?


texts:
I have two subfolders with Rasters named with the same parameters.
I am trying to figure out a way to Loop a multiplication of each raster into a folder with its pair of the other subfolder.
Looks like:.

SubFolder 1       SubFolder2         Result
178056_WGS.tif     178056_Clean.tif   178056_Mlt.tif
178057_WGS.tif     178057_Clean.tif   178056_Mlt.tif
178058_WGS.tif     178058_Clean.tif   178056_Mlt.tif
178059_WGS.tif     178059_Clean.tif   178056_Mlt.tif

Another option is to have all Raster in the same folder (workspace) but it is not at all desirable. Do you have any suggestion to do it in Arcgis 10.1 with Python (arcpy)?


texts: Selecting lines where start and end points within polygon using Virtual Layer

texts:
In QGIS 2.18 project there are several layers, lines and polygons respectively

My task is to select all lines when their start and end points are laying within polygons using Virtual Layer, see image below.

I can achieve the result that I want with the following query

WITH start_p AS (
SELECT l.id AS start_id
FROM line AS l, "Layer_2" AS p
WHERE ST_within(start_point(l.geometry),p.geometry)
),

end_p AS (
SELECT l.id AS end_id
FROM line AS l, "Layer_2" AS p
WHERE ST_within(end_point(l.geometry),p.geometry)
)

SELECT *
FROM line AS l, start_p, end_p
WHERE l.id = start_p.start_id
AND l.id = end_p.end_id

My Question: Why I can not achieve the same result using the following query

SELECT l.*
FROM line AS l, "Layer_2" as p
WHERE ST_within(end_point(l.geometry),p.geometry)
AND ST_within(start_point(l.geometry),p.geometry)

However, a query like (see below) will work

SELECT l.*
FROM line AS l, "Layer_2" as p
WHERE ST_within(start_point(l.geometry),p.geometry)


texts:
I would like to use a MAX SQL query for a column in a definition query, where there records that are being queried are the MAXIUM from the field UPDATE_DATE_DATE while grouping by the field NUMBERCYLA.  I am using this in an an SDE enabled geodatabase on SQL server 2008.

Here is what I have so far
[UPDATE_DATE_DATE] in (SELECT max( [UPDATE_DATE_DATE] ) FROM SO_SC_1 GROUP BY [NUMBERCYLA])


texts:
I have QGIS loaded on a couple of (Linux) computers, but by and large have used it as little more than an improved "Google Earth" with attempts at getting live data updates etc.; I've not seriously investigated scripting or advanced facilities.
I have two locations, one in the UK and the other in Germany. How can I plot on a map the line (a hyperbola?) that represents all points which are 100 km nearer the UK site than the German one?


texts: QGIS "transform error" messages - how to stop them?

texts:
In QGIS, I am getting an error message every time I edit a node. Also, there is a lot of lag in performance. The message is:

"Transform error caught: forward transform of
(10990.308086, 100347.255186)
PROJ.4: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +to +proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Error: latitude or longitude exceeded limits"

I hae been very careful to Set Layer CRS on all layers, to GDA94, or, in the case of the two raster layers, GDA94/MGA zone 54 (which is also the project CRS). This should all work. Shouldn't it?

But it is odd, that when I was using my legacy MapInfo TAB layers, they would show up in QGIS with weird projections that were not the projection in the TAB files, like some old Dutch projection and, more usually, NAD projections. I saved those layers as SHP files (resetting the projection to GDA94) to see if QGIS works better with SHP. And it did, for a few hours. Now the error messages are back. What is wrong? Is this a QGIS bug? Or is there a backoffice setting that I've missed?

EDIT: ADDITIONAL INFORMATION.

Below I provide the sequence of steps to create the transform error from a clean start.

Open QGIS
Set these options for Project CRS:

"Automatically enable OTF if layers have different CRS"
Do not "enable OTF by default"
Start new projects with GDA94
"When new layer is created, use project CRS."

Add raster layer 76949_zone54_mga.ecw (www.ga.gov.au/cedda/maps/1058)
Look at Project properties, witness that "OTF by default" is disabled, and CRS has changed to greyed-out MGA94 zone 54 (the CRS of the raster file)
Add shapefile layer (graticule from "Hamilton Special" GA 250K topo), downloaded from GA website.

Note that its .prj file contents are: GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
That is, it has a CRS of GDA94. The file loads and displays correctly. Look at project properties, and witness that it has automatically enabled OTF projection.

Create new shapefile point layer, with CRS selected as GDA94
Edit new shp layer, by making two points.
Exit edit mode and save.
Look at message log: two error messages logged with following info:

Transform error caught: forward transform of
(0.000000, 0.000000)
PROJ.4: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +to +proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Error: latitude or longitude exceeded limits

That is a replicable error, and all files used are pristine Geoscience Australia product, unquestionably ok. The fault lies somewhere in QGIS or in the settings used.

Yes I could create the new file in a CRS of MGA_zone54 (to follow the raster projection in this case), and that will work ok. But I would prefer that my vector files that aren't tied to a particular UTM zone. Anyway, GDA and MGA zone 54 are based on the same datum and ellipsoid, they are two sides of the same coin. MGA is just the projection of GDA94 longlat data.

So what is the problem?