texts: Openlayers : Zoom In and Map Layers Blocked wih WHITE space

texts:
I've a problem where something is blocking my map layer when I zoomed the map
click to view

function load() {
    map = new OpenLayers.Map({
        div: "map",
        maxExtent: new OpenLayers.Bounds(101.076797485352, 3.60602392425537,101.077194213867, 3.60696087346802)

        });
    // setup tiled layer
    floorplan7 = new OpenLayers.Layer.WMS(
                "assetmapper:floor_plan_t7 - Tiled", "http://192.168.56.101:8080/geoserver/assetmapper/wms",
                {
                    LAYERS: 'assetmapper:floor_plan_t7',
                    STYLES: '',
                    format: 'image/png',

                },
                {
                    buffer: 0,
                //    displayOutsideMaxExtent: true,
                    isBaseLayer: true,
                    yx : {'EPSG:4326' : true}
                }
        );

map.addLayers([floorplan7]);
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.addControl(new OpenLayers.Control.ScaleLine({}));
    map.zoomToMaxExtent();

}

What should I do to clear the 'WHITE' area?


texts: Google Earth Engine: Javascript API vs. Python API

texts:
Google Earth Engine seems to be most commonly used through either the Python API, or the Javascript API on the browser. Trying to decide for one, I have come across this comparison, from a Tyler Erickson presentation:

While some differences are clear (e.g. many more plotting options in Python), others are more abstract to me.

Questions:

In practical terms, what are the advantages of being able to make a "series of calls to Earth Engine" in the Python API?
Are there other differences not listed here, perhaps in regards to efficiency when exporting data, or the speed of computations?


texts: Getting AttributeError: 'DataFrame' object has no attribute 'to_file' from GeoPandas even though file converted into GeoDataframe instead of Dataframe

texts:
I have turned some R scripts for data cleaning/standardization of shapefiles into Python scripts.
Here is the R code:
x<-c(farm,field,"trialyield",2017)
file<-paste(x, collapse="_")
yield <- readOGR(".", file)
#identify column of dry yield
yield.df<-data.table::as.data.table(yield)
yield<-yield[,4]
names(yield@data)[1] <- "yield"

and here is what I did to try and convert the block of R code in Python:
y = farm + ' ' + field + ' ' + 'trialyield 2017'
y = y.replace(' ', '_')

trialyield = gpd.read_file(y)

title_col = trialyield.columns[0]
name_map = dict(zip(trialyield.columns[[4]], ['yield']))
trialyield.rename(columns=name_map, inplace=True)
trialyield[[title_col, 'yield']].to_file('trialyield_output.shp')

However, I get the current attribute error:
Traceback (most recent call last):
  File "./sommer_uofiaddington1_cleanandagg_2017.py", line 48, in <module>
    trialyield[[title_col, 'yield']].to_file('trialyield_output.shp')
  File "/opt/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py",
line 3614, in __getattr__
    return object.__getattribute__(self, name)
AttributeError: 'DataFrame' object has no attribute 'to_file'

if I recall gpd.read_file returns a GeodataFrame.
Why is it returning a DataFrame?


texts:
I have a set of LineStrings (streets) and a set of rectangles (atlas pages extents) that could overlap each-other.

And now I want to create a legend for an atlas like

street1, pages: 1,2
street2, pages: 2,3

or

page1, streets: 1,2,3,5
page2, streets: 2,4,5,6

Page numbers are rectangles' attribute. So I need to spatially join attributes, but relations between streets and pages are many-to-many and QGIS Spatial attributes join allows me only two options: write first occurrence or aggregate over attributes values.

In other words, I want to get the analog of Postgis spatial join like below and write it out (in .csv or other format):

select page.number, street.number, street.name
from
   street, page
where
   ST_INTERSECTS(street.geom, page.geom)

I could run a python script that will iterate over pages, take streets, and write cross joined attributes into a .csv file. But maybe somebody already done this before and has published a solution?


texts: QGIS joined attributes not visible in categorized symbology

texts:
I am trying to make a thematic layer using the attributes from my now main table which were previously joined.

I would like it to look something like this:

As visible in picture, the columns starting with Gotenburg, Gothenbu_1 etc are in the table joined.

But when I try to create categories based on its values, I do not have the option to use those attributes. As seen in next photo.


texts: I am using QGIS 2.8.3.After dissolving a shapefile i will check for geometry validity there it shows many green cross errors( two nodes overlap at a time which results in that green crosses) i got some 200 crosses like this. How do i can remove this automatically because doing it manually consume lot of man hours


texts: Assigning intersecting line length dependent value to raster during polyline conversion using ArcMap?

texts:
I am converting a road network polyline into raster format- my polyline shapefile has an exhaust emissions rate attribute which is in kg/km

What I want to achieve is a raster that has (at say 5m resolution) cells with computed emission values that are a function of the length of road that falls within the cell times the emission rate.

For example, for a cell of 5m if the road line intersects at 45 degrees through each corner of the cell the raster value should be the emission rate times 7.1 (a^2+b^2=c^2)- at the moment ArcMap is just taking the raw emission rate and making no correction for geometry.

Can anyone help?


texts:
I'm working in a relatively large project in QGIS 2.14, mostly loaded with .shp files. The issue is that lately got painfully slow when zooming or changing map view. I've tried to delete the largest layers and performed spatial indexing for all layers with no luck.

Is there a way to identify which layer is using more resources and slowing down the project?


texts: Speeding-up processing of pair-wise distance calculation in Python?

texts:
I am trying to do a cluster analysis using DBSCAN for my time series NDVI image in Python. I am using distance time warping (DTW) to measure distances between my time series. With this, I am computing pairwise distances using DTW which will be eventually be an input to DBSCAN. My problem, pairwise calculation seems very slow. I browsed some techniques, such as multithreading, GPU programming, splitting data into blocks. I would like to ask, which could be the approach to speed-up my processing for calculating pairwise distances. An approach that I can simply implement on a 2-core x 8GBRAM x i7 processor laptop.

Some information:

The input numpy table (mod_hypercube.npy) is an output of image stacking which I did in my other script. It is an image stack stored as numpy table.

The DBSCAN parameters are still experimental, will explore different combinations when I have the pairwise metric.

The numpy table dimension is

In [4]: t.shape
Out[4]: (138L, 624L, 645L)

The code is:

import numpy as np
import math,gc
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN

def DTWDistance(s1, s2):
    DTW={}
    for i in range(len(s1)):
        DTW[(i, -1)] = float('inf')
    for i in range(len(s2)):
        DTW[(-1, i)] = float('inf')
    DTW[(-1, -1)] = 0
    for i in range(len(s1)):
        for j in range(len(s2)):
            dist= (s1[i]-s2[j])**2
            DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], DTW[(i-1, j-1)])
    return math.sqrt(DTW[len(s1)-1, len(s2)-1])

t = np.load("E:\\file\\path\\mod_hypercube.npy")

distMatrix = np.empty(shape=(t.shape[1] * t.shape[2],t.shape[1] * t.shape[2]))

# this is the dimension of the hypercube matrix
idx = 0
for r in range(0,t.shape[1]):
    for c in range(0,t.shape[2]):

        X = t[:,r,c]
        idx_2 = 0
        distances = np.empty(shape=(t.shape[1] * t.shape[2]))
        for j in range(0,t.shape[1]):
            for k in range(0,t.shape[2]):
                y = t[:,j,k]
                distances[idx_2] = np.round(DTWDistance(X,y),2)
                idx_2 += 1

        distMatrix[idx] = distances
        print "this is idx1=%s" %idx
        idx += 1

np.savetxt("E:\\file\\path\\distMatrix.csv", distMatrix, delimiter=",")
db = DBSCAN(eps=0.3, min_samples=10,metric="precomputed").fit(distMatrix)
# parameters=‘auto’, ‘ball_tree’, ‘kd_tree’, ‘brute’
#(eps=0.5, min_samples=5, metric=’euclidean’, metric_params=None, algorithm=’auto’, leaf_size=30, p=None, n_jobs=1)
print db.labels_

plt.imshow(distMatrix)
plt.show()


texts: Why are all tiles white when using the OpenLayers plugin?

texts:
I installed the OpenLayers-Plugin and set the CRS to EPSG:3857.

When I go to Web / OpenLayers plugin I get the menu with the various maps, but whatever map I choose, all the tiles are white.
The layer shows up in the legend correctly, if there is another layer underneath, it is hidden by the white tiles.

When I open the OpenLayers toolbar, the small preview map is also white.

What could be the reason?

edit:
The shapefile I used can be found at
www.biooffice.at/java/adscita_geryon_italy.zip

I am adding screenshots with OpenLayers active and off:


texts: Editing only selected features in ArcMap?

texts:
I've got a big features layer with more than 30 000 points. I want to edit some of them using the snap tool, but I'm afraid it's going to edit some points I don't want to be edited.

Is it possible to only edit selected features?


texts: Create new layer by selection within ModelBuilder in ArcMap 10?

texts:
I have a basic abutters model built so far. Right now I have it set up to:

Select a parcel by an expression.
Buffer parcel by a value.
Select by location which highlights all parcels within the buffered area.

I would like to have a new feature class created with only the selected parcel between steps 1 and 2. I tried putting copy features between the output of the first selection and the buffer, but a new feature class of just the selected isn't created.

I would like to do this so the selected parcel has a separate symbology than the rest of the parcels.


texts: What are some intuitive ways to have fuel optimal routes?

texts:
Given a road network and two points: start and end, what are some good reliable? ways to estimate the fuel consumption over a road segment (i.e. a network edge).

My goal is to output the fuel-optimal path ( the one with least fuel consumption among all possible paths) instead of shortest/fastest path. I am working on the OSM raw data only and I have no access to 3rd party APIs/software.

I found it a bit tricky specially when the there is no information over signals and turns. All I have is the vehicle info and the OSM data.

A naive way: the longer the segment the better with carefully chosen penalty so we don't end up with the longest path. Possibly also by incorporating the straightness of the road but don't know how to estimate it given its start and end coordinates. Any suggestions on estimating fuel-consumption for a road segment?


texts: Displaying landcover classes using Google Earth Engine

texts:
I was wondering if theres a way to display single Landcover classes from the Corine Land Cover-Dataset in GEE. I know i somehow have to select the landcover_class_names or values in the property section of the image but I'm not sure how to do that.

var corine = ee.ImageCollection('COPERNICUS/CORINE/V18_5_1/100m').first();

print(corine);
print(corine.propertyNames());
var lc_value = corine.get('landcover_class_values');
print(lc_value);

Thats my code so far. Now i only want to display the forest areas.


texts:
I have a KML from Indy Neighborhoods ArcGIS. I have the results of a factor analysis based on neighborhood traits from savi.org. From this, I did a clustering. I would like to combine the contents of the clustering outcome spreadsheet, consisting of cluster, color, and scores of factors, to the KML.
How would I do this?


texts:
I'm trying to create a 8 bit GeoTIFF with this code.

Is there a way to use unsigned 8 bit? I only found the GDT_UInt16?

import numpy as np
from osgeo import ogr, osr, gdal
from osgeo.gdalconst import *

dataset = gdal.Open('test.tif')
driver = dataset.GetDriver()
driver.Create("out.tiff", cols, rows, 1, GDT_UInt16) # <- here

out_band = out_dataset.GetRasterBand(1)
out_data = np.zeros((rows,cols), np.int16)
out_band.WriteArray(out_data, 0, 0)
out_band.FlushCache()
out_dataset.SetGeoTransform(dataset.GetGeoTransform())
out_dataset.SetProjection(dataset.GetProjection())


texts: Problem with an Expression in Attribute Assistant?

texts:
I am using ArcGIS 10.3.1 for Desktop in a versioned environment.  I am also using ESRI Attribute Assistant to autopopulate fields using a custom expression.

My Table Names, Field Names, Value Method, etc... all are working fine.  My problem is with my expression below:

iif(InStr(UCase([REQUEST_TYPE]), "SEW") > 0,  VbNull , "Inspection")

The problem here is the expression, if it finds "SEW" in the [REQUEST_TYPE] it is returning a 1 instead of a <NULL> , otherwise it is returning "Inspection" fine.

Why is it returning a 1 ?  Also how can I fix this.  I have tried 'Nothing' , Nothing , and <NULL>  instead of VbNull but nothing seems to work.  What would the correct Expression be?

I tried some more things,  I currently have the EXPRESSION setup as:

iif(InStr(UCase([REQUEST_TYPE]), "SEW") > 0,  "" , "Inspection")

so that it returning a "" if "SEW" is found in [REQUEST_TYPE] , though this is better, I would still prefer a NULL value if possible.


texts: Making custom markers resizable in Leaflet?

texts:
This fiddle (from Stefan_Fairphone) is a quick example of what I need to do for a project. But i need to add two features for these custom markers: make them scalable and spinnable.

I'm trying to use the L.Path.Transform plug-in, but it doesn't seem to work, I always get the error L.marker(...).transform.enable(...) is not a function OR Can't read the property of 'enable' of undefined

I also had a script Here that uses the ID of clone, but it doesn't work with what I have now (don't mind the infinite cloning), the error I have is L.marker(...).attr(...) is not a function

And if I use Leaflet.StylEditor, the custom marker will turn into a regular marker.

I'm trying to make my marker resizable by modifying the iconSize value by doing this :

target = document.getElementById("mapid");
        target.ondragover = function (e) {
            e.preventDefault();
            e.dataTransfer.dropEffect = "move";
        };

        target.ondrop = function (e){
            e.preventDefault();
            imagePath = e.dataTransfer.getData("text/plain");
            coordinates = mymap.containerPointToLatLng(L.point([e.clientX,e.clientY]));
            function resize(){
                $(document).keydown(function(event){
                    width = 60;
                    height = 40;
                    Size = [width, height];
                    if(event.which == 83){ //S key
                        height+=5;
                    };
                    if(event.which == 90){ //Z key
                        height-=5;
                    };
                    if(event.which == 91){ //Q key
                        width-=5;
                    };
                    if(event.which == 68){ //D key
                        width+=5;
                    };
                    Size = [width, height];
                    return Size;
                });
            };
            Size = resize();
            L.marker(coordinates, {icon: L.icon({iconUrl: imagePath, iconSize: Size}),draggable: true}).addTo(mymap);
        };

But it's not working...

I'm looking for a way/workaround to make custom markers, resizable and make them spin around.

After a couples of tries, I realized I could use th attr function a little different

Made a good step forward : I'm trying to combine both of the scripts, since they complete each other for what I need, and using a class element leaflet-marker-icon with the attr, I can use it to set the ID of my marker, like so :

$('.leaflet-marker-icon').attr('id', 'rotate');
var angle = 0;
var width = 60;
var height = 40;

$(document).keydown(function(e){
    if(e.which == 65){ //"A" key
        angle-=5;
        left();
    };
    if(e.which == 69){ //"E" key
        angle+=5;
        right();
    };
    if(e.which == 83){ //"S" key
        height+=5;
        smaller();
    };
    if(e.which == 90){ //"Z" key
        height-=5;
        bigger();
    };
    if(e.which == 81){ //"Q" key
        width-=5;
        thinner();
    };
    if(e.which == 68){ //"D" key
        width+=5;
        wider();
    };
});

function left(){document.getElementById('rotate').style.transform = "rotate("+angle+"deg)";};
function right(){document.getElementById('rotate').style.transform = "rotate("+angle+"deg)";};
function smaller(){document.getElementById('rotate').style.height =  height + 'px'};
function bigger(){document.getElementById('rotate').style.height = height + 'px'};
function wider(){document.getElementById('rotate').style.width = width + 'px'};
function thinner(){document.getElementById('rotate').style.width = width + 'px'};

So, the resize works fine, but the rotation is kinda weird. during the rotation sequence, the marker lose the translate3d and its coordinates, and goes to the to last dropped place, and when an action occured (zoom, unzoom, drag), the marker lose its angle value and goes back to its original area, no idea why

And, I'll need to swap that ID attribute, cuz I'll have a good dozens of custom markers inside the map...


texts: Can GeoServer Mosaic extension footprints be something other than a shapefile?

texts:
According to the GeoServer documentation for the Image Mosaic extension...

http://docs.geoserver.org/2.4.x/en/user/tutorials/imagemosaic_footprint/imagemosaic_footprint.html

It asks for a shapefile as the footprint data set. Does anyone know if the footprint can be  any other type of data set, preferably PostGIS/PostgreSQL or ArcSDE / SQL ?


texts:
I have a selection of python plugins (one rather large one) that I have developed under Qgis 1.8 and that I am migrating to qgis 2.0.  The error that the Plugin Manager gives me when I load these plugins (This plugin was developed under 1.x) is rather unhelpful.  It also arises in situations that have nothing to do with the API Version.  For instance, a Network error yesterday introduced a stream of random characters into a few of my code files (Fun!) and this led to the Plugin Manager telling me my plugin was old. Eventually, I found out that a tab and a dollar sign had been inserted in front of a comment in my metadata.txt, and that was the problem.

Is there any way to see a more detailed error when QGIS 2 rejects a plugin as 1.x?

If not a proper trace, then at least what class, or even what file was rejected?

The instructions regards API changes at the Qgis Hub Wiki are not as useful as a simple error message.


texts: ArcGIS - Is It Possible To "Lock" Metadata

texts:
I've been asked by a customer if I can "lock" the metadata in their various geodatasets.  I immediately asked why, and was told that they don't want other organizations/users editing/deleting certain tags.

A quick search hasn't really turned anything up.  Anyone have any ideas/suggestions on how I could accomplish this?


texts: I have to convert a coordinate from UTM to MGRS. I think I have the most of it covered, but I have some difficulties to calculate the correct Latitude-band. How can I decide which band the coordinate has?


texts: How to get TOC and Widgets (scalebar, overviewmap and basemapgallery) all working with ArcGIS JavaScript API?

texts:
Yesterday I spent many hours trying to get a scalebar, overviewmap and basemapgallery working in one webmap using nilu'S TOC.

Unfortunately without success :(

I have one version, where scalebar and overviewmap work (but not the basemapgallery) and another version, where scalebar and basemapgallery work (but not the overviewmap).

Here is a fiddle of the version with working scalebar and basemapgallery work (but not the overviewmap):

As you can see in the fiddle the OverviewMap Widget is already included.

I have absolutely no idea what I have missed in the code. However, I think something must be misplaced.

Can anyone suggest what may be astray?


texts: Georeferencing historic airphotos along constantly-moving section of San Andreas fault?

texts:
I'm a geology graduate student working on a master's thesis project in which I'm measuring creep (movement) over time along a section of the San Andreas fault (in California, USA) using SfM derived from overlapping historic aerial photos. The airphoto series' I'm using are from various time periods spanning the last ~80 years and none of them are georeferenced. The section of the fault I'm focused on is the creeping section, which means it is constantly slipping - at an average rate of about 30 mm/yr. Because both sides of the fault are always moving (in opposite directions relative to one another), any potential ground control points from one time period will not match spatially with GCPs from any other time period.
The end goal of the project is to compare 3D point clouds (from SfM, using Agisoft Metashape Pro) from each time period to one another for eventual ICP analysis to measure the respective offsets. The only way this can work is if the data from each time period are precisely georeferenced to one another. Because I have no stationary GCPs over time, I am at a loss as to how to accomplish this. Since the data does not necessarily have to be georeferenced to the Earth (although that would be ideal) and only have to be georeferenced to one another, one idea I had was to use GCPs on only one side of the fault; but that seems likely to produce unacceptably high levels of error. I'm also wondering if it may be easier to just produce the point clouds as-is (without any coordinate system) and deal with georeferencing those rather than the airphotos.


texts: Converting shapefile in NAD27 from meters to decimal degrees to match xy data?

texts:
I have a shapefile of a boundary (polygon) that has the following projection:

  NAD_1927_UTM_Zone_11N
  Projection:   Transverse_Mercator
  False_Easting:    500000.00000000
  False_Northing:   0.00000000
  Central_Meridian: -117.00000000
  Scale_Factor: 0.99960000
  Latitude_Of_Origin:   0.00000000
  Linear Unit:  Meter

  Geographic Coordinate System: GCS_North_American_1927
  Datum:    D_North_American_1927
  Prime Meridian:   Greenwich
  Angular Unit:     Degree

I need to pull climate data (in netcdf files) that is within this boundary.  That data is in x y coordinates and has the following information:

  XLAT
             Size:       219x285
             Dimensions: londim,latdim
             Datatype:   single
             Attributes:
                         units = 'deg'
      XLONG
             Size:       219x285
             Dimensions: londim,latdim
             Datatype:   single
             Attributes:
                         units = 'deg'

There's no other metadata, so in order to align these files, I need to convert the shapefile to decimal degrees.  I tried using the Project tool, selecting the same coordinate system, but I don't see decimal degrees in the drop-down.


texts: Issues ST_MakePoint in PostGIS with Geoserver

texts:
I am trying to setup a GIS system to use PostGIS and Geoserver. I have loaded all of the data into a table. The data that is received has long/lat coordinates in Degrees, Minutes, Seconds, Direction, each of which is a separate column. What I did is created a column (location) of type geography and I'm trying to give proper info for each record. The command that I am using is this (of course this is only in the Western Hemisphere):

UPDATE *tablename* SET location = ST_SetSRID(ST_MakePoint((lat_degrees+(lat_minutes+lat_seconds/60)/60), (-1*(long_degrees+(long_minutes+long_seconds/60)/60))), 4326);

My DMS -> DD conversions should be correct, and I've verified that by creating separate columns to test each individually. But the problem is, my geoserver doesn't detect the proper bounds. When I manually set the bounds (I set them to -180,-180,180,180) and check my data in Google Earth, it only zooms to 0,0 and I can find some horribly warped data, in Antarctica.

One thing to note... not all of my records include location data, could that be an issue? Or is ST_MakePoint simply not creating the right location for my data?

Thanks for the help!


texts: How to get the data from Google Earth Engine?

texts:
I want to initialize three variables: one for water, another for buildings and the last one for the other things for 2010,2015 and 2020 of a certain location starting from the month of March using Landsat 8 Tier 1 Collection 1 OLI Raw scenes. How can I do this?
This much I have done and being a beginner to this field I don't have much experience working with this.


texts: How to add map layer of Italy in map canvas in standalone QGIS application

texts: I have map canvas and I need to add map layer of Italy. I already succeed to add map of whole world, but when I zoom it to Italy I get nothing. I do research and I found this page: http://mobac.sourceforge.net/forum/viewtopic.php?t=519 with XML code. But I am not sure how to add that. When I added map of whole world I just put URL in QgsRasterLayer and that was it. I find i cookbok in cookbook some guidlines, but nothing specific.


texts: There is combine function is ArcGIS which combine multiple raster value attributes in one raster layer. I would like to know what's the same process in QGIS?


texts: ArcMap to automatically create polygons from shapes that are present in an image raster

texts:
I have a geolocated greyscale geojpeg image showing a network of geological fractures in the ground. They show up as bright lineations and create a reticulate pattern (think giraffe pattern). I've converted the image into a binary where the raster values of the fractures are 0 (black) and the remaining pixels not associated with the fracture pattern are white. The fractures are between 1 and 5 pixels wide. There are also some noisy pixels not associated with the fractures.

How do I convert this binary image into a vector file? I want to end up with a polyline shapefile of the fracture network so I can quantify some metrics, like the angle the fractures make to each other, etc. This would replace having to digitize the pattern by hand.

I've tried the raster to polyline and polygon tool but it creates vectors around each individual raster pixel, which is not what I want. I have the Spatial Analyst and 3D Analyst extensions.

Any tips on tools or workflow?

Update: screenshots attached


texts: Raster on-the-fly projection in ArcGIS has different resolutions

texts:
I have a map document authored in ArcMap.  The spatial reference for each of the raster layers is Alaska State Plane.  If I set the Data Frame projection to Web Mercator (Auxiliary Sphere), then my raster layers have a resolution of about 2 meters, when the actual data has a resolution of about 2 feet.  I'm trying to publish the map document as a Map Service (ArcGIS Server 10.0) but the quality of the imagery is not acceptable for the use case.

Is there any way to make the imagery draw at 2-ft. resolution, or do I need to reproject the data (about 500 GB of raster imagery)?


texts: Cross-validation of TPS and kriging outputs from the fields package in R

texts:
I am looking for existing package or functions that will allow me to cross-validate multiple outputs of Thin plate spline (TPS) and Kriging from the fields package.

Does anyone know of any available package or function in R that would do this?


texts: 1944 Cyprus Grid coordinate data projection problem in MapInfo

texts:
I have coordinate data from the 1944 British Army General Staff, Geographical
Section (GSGS) Cyprus map which was cast on the Cassini projection, but was referenced to the Clarke 1880 ellipsoid where a = 6,378,249.145 m and 1/f = 293.465.

I gather from Clifford Mugnier's 2006 article that the latitude of origin should be (φo) = 35º 00’ 00”N and the central meridian (λo) = 33º 19’ 00”E, with a False Easting = 200 km, and a False Northing = 150 km.

I've added these data to the prj file for the Cassini-Soldner projection ("Cassini-Soldner Cyprus", 30, 1000, 7, 33.316667, 35.00000, 200000, 150000), which puts the points in roughly the right location, but their locations are distorted slightly, resulting in some sites appearing in the sea, etc.

Is there any way of getting the points in the exact right location? Have I forgotten to include some information?

Thanks


texts: Showing contour line elevation labels in QGIS?

texts:
I have a contour map in DXF imported to QGIS. How do I turn on the elevation labels for each contour so they are visible. I thought I had done it by clicking on the layer and selecting properties > labels - single label > Label with - layer. All I get is a '0' for each contour.

Is there info in the DXF missing?


texts:
Anyone knows how can i open a .rrd file (should be a raster)in QGIS 2.6? Comes with more 3 files.
Should be any difference using it in OS or Windows?


texts: Removing features whose labels are not displayed?

texts:
Using ArcMap 10.0, I am working on a map that shows many villages. Each village has a name attribute that I want to use for a label, but in some areas there are so many villages that not all can be labeled. I am using the Maplex label engine and have not checked the "Never Remove (allow overlap)" box.

I am quite happy that not all village names are displayed as this would be way too much. And it is ok that Maplex chooses which labels to display, as I have no way to know which village is more important than another one. This also means that I can't base the labeling on a field such as population or something else, even though that would be less random.

However, I would also like the village points that are not labeled to similarly not be displayed. I think it unlikely that someone would want to see the point without knowing the name of the village. If all points are shown, it is also difficult to tell which point a name label goes with in crowded areas. So instead of choosing just labels not to be shown, I also want the computer to chose features not to be shown.

To give you a better idea on what I mean:

Is there a way to make these features/points not draw if their label is not shown?


texts: How to add CSS styling in Mapbox GL Popup?

texts:
I am aware this is similar to many other questions, but I have been unable to determine how to add the correct code to style my popups in mapbox-gl. I'm very new to CSS, Mapbox-gl, and coding in general. I would like to be able to add banners, change the colors and fonts of different sections, and other additions. I've read this and several posts on SE, which pointed me in the right direction, however I have still been unable to make this work.

Here is a sample of my code - I think I know where I am doing something wrong but I am not sure what to do right.

<meta charset='utf-8' />
<title>Points on a map</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.31.0/mapbox-gl.css' rel='stylesheet' />


texts: When creating polygon features in ArcView 10, using the auto-complete polygon tool, I create a feature neighboring another polygon and ArcView hangs or processes for over 15 seconds. The feature is successfully created, but it takes forever. The polygon layer being edited is in a file geodatabase (although location of the feature, in or out of a gdb, being edited appears to have no bearing on the problem). I've done the Normal.mxt rename and all the tried-and-true solutions and ESRI has been of little help.


texts: Merging polygon slivers to largest adjacent polygon using ArcGIS Desktop?

texts:
I've got several hundred thousand slivers in my dataset as a result of cutting in a layer.

Does anyone have a good solution as to how to have ArcGIS Desktop evaluate the position of each sliver and determine where they need to be allocated?

I want each sliver to be assigned to the largest adjacent polygon that it shares more than one node with.


texts:
I have a shapefile where the geometry is a mess, lots of gaps, slivers, and holes. I need to correct these geometry errors, but I am having trouble with the process.

I have read the pages here, here and here, and looked at the SE questions:

Repairing geometry using QGIS?,
Cleaning polygon after using dissolve to repair geometry in QGIS? and
Fixing geometry validity errors in QGIS?

amongst others, and the GRASSwiki pages.

However, I am still uncertain about the steps I should be taking to correct these, particularly the best workflow and the settings I should use.

I can't seem to get anywhere with Dissolve or Eliminate Sliver Polygons, and so I have been experimenting with GRASS in QGIS.

Currently, my steps are;

Check for errors on the file using Vector > Geometry Tools > Check Geometries, (I have tried with Topology Checker but it keeps crashing Q). I have been using these options:

(Check Geometries says there are 3330 errors (in a shapefile with only 1900 polygons). I have tried fixing some of these in the Results tab of Check Geometries, but it is very slow and each time I resolve one error it finds 1000+ more. So I have been attempting to use GRASS instead)

Steps I'm using

Open GRASS mapset and import the shapefile using v.in.ogr
I opted for snapping of 9 because that seems to preserve the data. The ‘pixels’ are 10m x 10m, if I increase the snapping tolerance to 10 some of the gaps begin to fill in but the square pixels begin to warp, and some of the polygons are no longer imported.
Basic import into GRASS and export again (with minimum area of 101) leaves 2330 errors. So start trying to improve the errors using v.clean  (should I untick centroids for any of these?) bpol, break, rmdupl, rmarea (threshold = 101, ie. > one pixel).
Export with v.out.ogr (with 'centroids' unticked).
Check the shapefile with Check Geometries again. Only removing the smallest polygons, (either during import or with 'rmarea') seems to have any effect.

I’m left with 2325 errors, all either ‘Polygon with hole’ or ‘Sliver polygon’. This still seems like a huge number of errors. The ‘sliver’ polygons are small, but they are not 'errors' in the sense that I want to retain them for future analysis. And can polygons really not have holes? I don’t really want to have to delete any more than I absolutely have to. The shapefile is based on a landcover classification, and if I start to alter the image too much it will lose its value.

So, to sum up:

Are these steps the right ones?
If they are, why have I still got so many errors?
And, how can I resolve them without losing too much information?


texts:
I am currently developing a QGis plugin and I need to create and remove directories and I use PyQt.
I don't know if my problem comes from PyQt or QGis. So my program creates folders, that works all good smoothly and stuff wow cool XD, then I can remove the empty folders with rmdir from QDir, that works too. the problem comes when there is files in the folder... In order to remove the non empty folders I first remove the files within using the QDir's remove method, that works fine all the files burns in hell and when I want to remove the folder, the rmdir method return True, like if everything was perfect but in fact no XD the folder's still here and became like "blocked" and I can't remove it from the program nor from windows "manually" ! And then only when I exit QGis, all the deletions are applied... I don't know where it comes from... I guess it's because of something like the directory is used by QGis but I don't know... so if you want to see the code ask for it because it's big and I didn't want to throw evrything ^^


texts: Updating and merging geometries using a 3rd table as the connector

texts:
I'm looking for some help with the following.
I have 3 tables in my Oracle Database.
Table A:
IDnumber  Geometry
001       Geom1
002       Geom2
003       Geom3

Table B:
IDnumber  Fieldnumber (neither of these are unique fields)
001       Field1
001       Field2
002       Field3
003       Field1
003       Field4

Table C:
Fieldnumber  Geometry
Field1       Geometry1
Field2       Geometry2
Field3       Geometry3
Field4       Geometry4

I want to update the geometries from Table A with the geometries from Table C using Table B as the connector. If there is more than 1 Fieldnumber for each IDnumber in Table B, I want the geometries from Table C to merge.
My update SQL so far:
Update TableA a
Set
a.Geometry = (select c.Geometry from TableC c where c.Fieldnumber = b.Fieldnumber)
Where
a.IDnumber in (select b.IDnumber from TableB b)

This only works when the TableB IDnumbers are unique. How can I merge them? With StUnion? If so how does that work? I don't understand the code.


texts: Why am I getting a syntax error at "RETURN result;" when creating a CartoDB pl/pgsql function?

texts:
Despite CartoDB appearing to run PostGIS 2.2 (see below), it still doesn't support multi-point geometries in ST_Split.

SELECT postgis_full_version();

POSTGIS="2.2.0" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012"
GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.7.8"
LIBJSON="UNKNOWN" (core procs from "2.2.0" need upgrade)
RASTER (raster procs from "2.2.0" need upgrade)

So I tried to use the workaround function from this answer. However, though creating that function works on my local machine, I get the following error when I create the function in cartoDB:

syntax error at or near "result"``

Here's the function I'm creating:

DROP FUNCTION IF EXISTS split_line_multipoint(input_geom geometry, blade geometry);
CREATE FUNCTION split_line_multipoint(input_geom geometry, blade geometry)
  RETURNS geometry AS
$BODY$
    -- this function is a wrapper around the function ST_Split
    -- to allow splitting multilines with multipoints
    --
    DECLARE
        result geometry;
        simple_blade geometry;
        blade_geometry_type text := GeometryType(blade);
        geom_geometry_type text := GeometryType(input_geom);
    BEGIN
        IF blade_geometry_type NOT ILIKE 'MULTI%' THEN
            RETURN ST_Split(input_geom, blade);
        ELSIF blade_geometry_type NOT ILIKE '%POINT' THEN
            RAISE NOTICE 'Need a Point/MultiPoint blade';
            RETURN NULL;
        END IF;

        IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN
            RAISE NOTICE 'Need a LineString/MultiLineString input_geom';
            RETURN NULL;
        END IF;

        result := input_geom;
        -- Loop on all the points in the blade
        FOR simple_blade IN SELECT (ST_Dump(ST_CollectionExtract(blade, 1))).geom
        LOOP
            -- keep splitting the previous result
            result := ST_CollectionExtract(ST_Split(result, simple_blade), 2);
        END LOOP;
        RETURN result;
    END;
$BODY$
LANGUAGE plpgsql IMMUTABLE;


texts: Tutorial on how to model flooding using data from gauging stations along rivers in QGIS?

texts: Can someone please direct me to a tutorial where I can learn more about how to model flooding using data from gauging stations along rivers using QGIS?


texts: Open and plot a large raster with Python xarray

texts:
I have an ASCII Esri file of the type:
ncols 5201
nrows 4001
xllcorner -65
yllcorner -35
cellsize 0.025
NODATA_value -9999
-9999 -9999 -9999 255 256 -9999 -9999 ...
-9999 -9999 -9999 -9999 -9999 -9999 -9999 ...
-9999 -9999 -9999 -9999 255 -9999 230 ...
...

I'm trying to import it with rioxarray by doing:
import xarray as xr
import rioxarray

da = rioxarray.open_rasterio('data.asc',masked=True)
da.plot()

What I get is a graph with the correct coordinates, but no data shown:

Do you know why this is?
da is imported as an xarray.DataArray and looks like this:


texts: Cannot clip because "unable to open database file". Do I need to create a geopackage from shapefiles, and if so, how?

texts:
I am working with some shapefiles in QGIS 3.2.1, and I would like to clip one shape file to another, and I am having the following error message.

Input parameters:
{ 'INPUT' : '.../Data/landcover/zwe_gc_adg.shp', 'OUTPUT' : 'test', 'OVERLAY' : '.../Data/cyclone_affected_areas/cyclone_affected_districts.shp' }

Could not create layer clipped.gpkg: Creation of data source failed (OGR error: sqlite3_open(clipped.gpkg) failed: unable to open database file)

I understand the problem is that I am using shapefiles, and I think I should be using a geopackage, which has layers, but I'm unsure what I need to do (I may just need a good introduction to creating a geopackage from existing shapefiles). Can I then work with the layers of the shapefile?


texts: Using round function in QGIS expression?

texts:
I created a grid of rectangular polygons with the Create Grid function, from the Vector Creation group of the Processing Toolbox in QGIS 3.4.2.
The rectangles measure half a degree of longitude by half a degree of latitude.

When I want to label the polygons, I write an expression that rounds the values of latitude and longitude attributes to 2 decimal places, since the attributes created in the grid brings many decimal places.

But the round function rounds the values badly.

I know that I can create new attributes that express the coordinates with fewer decimal places and label the entities based on those attributes, but I was struck by the behavior of the function. I expected it to round to the nearest value.

Does anyone know what the problem is that the round function wrongly rounds off the decimal number?

Attached screenshots of Source Fields and the Attribute Table.


texts: Using arcpy dictionaries to select attributes from a different layer

texts:
I'd like to incorporate the use of dictionaries to access data between tables. So far I have been able to create the dictionaries that I need from my feature class fields.  I have one dictionary that lists multiple Property values for one key [(u'1010305', [u'1014414', u'1014395', u'1014381']), (u'1010306', [u'1014396', u'1014415', u'1014432', u'1014414', u'1014395', u'1014382', u'1014381']).... I would like to use the values associated to each key to select rows ("select by attribute") in another table.  While I can call the values by key or item in the command window I can't seem to incorporate them into a SelectLayerByAttribute function.  This is the code so far :

# Purpose: Create a dictionary of values that combines Recipient Claims to Allocating Claims from the BufferRecipientClaims feature class.
# Usage: No arguments needed.
# Source : Python for ArcGIS, Laura Tateosian, 2015, Springer, Chapter 18 p.349.

import arcpy, numpy

arcpy.env.workspace = "CURRENT"
fc = "RecipientBufferClaims"

# Populate the dictionary,
# accumulate a list of Recipient Claims (Value) for each Allocation Claim (Key) type.
dict = {}
SearchCurs = arcpy.da.SearchCursor("RecipientBufferClaims",["RecClmNum","AllocClmNum"])
for row in SearchCurs:
    Recipient = row[0]
    Allocation = row[1]
    if dict.has_key(Recipient):
        dict[Recipient].append(Allocation)
    else:
        dict[Recipient] = [Allocation]
del SearchCurs

All of which works.  I'd like to use something like the following to select rows in a new table using the returned values for the key "1010305" but I'm getting an error message

for k, v in dict.values["1010305"]:
    AllocationVariable = AllocationVariable(v)
    print "the allocation claim numbers are  {0} ".format(AllocationVariable)
    arcpy.selectLayerByAttribute(fc,""" "AllocClmNum" = AllocationVariable """)

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ValueError: too many values to unpack

Edit: I did manage to get this to update another table, unfortunately it did so without filtering the values associated to the key 1010305, but I'm getting closer (code below).

with arcpy.da.UpdateCursor("ActualTbleCopy2SortedAllocationClaims",["AllocClmNum","CorrAllocCreVal"]) as UpdateCurs:
    for row in UpdateCurs:
        for k, v in dict.items():
            RecipientNumber = str(k)
            if RecipientNumber == ("1010305"):
                row[1] = 9999
                UpdateCurs.updateRow(row)


texts:
I'm trying to build a debug version of QGIS on windows env for this I'm using the bat file package-nightly, I followed all steps presents in INSTALL file but without success.

When I executed the bat file I received this error:

-- Build files have been written to: C:/Users/renan.b/QGIS/ms-windows/osgeo4w/build-qgis-dev-x86
CLEAN: 29/06/2017 14:36:16,58
Cleaning... 10 files.
ALL_BUILD: 29/06/2017 14:36:17,13
[1/4225] Building CXX object src\native\CMakeFiles\qgis_native.dir\qgsnative.cpp.obj
FAILED: src/native/CMakeFiles/qgis_native.dir/qgsnative.cpp.obj
C:\PROGRA~2\MICROS~1.0\VC\bin\cl.exe   /nologo /TP -DENABLE_TESTS -DNOMINMAX -DQGISDEBUG=1 -DQT_DISABLE_DEPRECATED_BEFORE=0 -DQT_NO_CAST_TO_ASCII -DWI
TH_BINDINGS -DWITH_QTWEBKIT -D_CRT_NONSTDC_NO_WARNINGS -D_CRT_SECURE_NO_WARNINGS -D_USE_MATH_DEFINES -Dqgis_native_EXPORTS -I. -I..\..\..\src\native -
Isrc\native /DWIN32 /D_WINDOWS /W3 /GR /EHsc -DSPATIALITE_VERSION_GE_4_0_0 -DSPATIALITE_VERSION_G_4_1_1 -DSPATIALITE_HAS_INIT_EX  /wd4100  /wd4127  /w
d4190  /wd4231  /wd4244  /wd4251  /wd4275  /wd4505  /wd4510  /wd4512  /wd4610  /wd4706  /wd4800  /MD /ZI /MP /Od /D NDEBUG /D QGISDEBUG -DQGISDEBUG=1
/showIncludes /Fosrc\native\CMakeFiles\qgis_native.dir\qgsnative.cpp.obj /Fdsrc\native\CMakeFiles\qgis_native.dir\ -c C:\Users\renan.b\QGIS\src\native
\qgsnative.cpp
c:\users\renan.b\qgis\src\native\qgsnative.h(34) : error C2065: 'default' : undeclared identifier
c:\users\renan.b\qgis\src\native\qgsnative.h(34) : error C2258: illegal pure syntax, must be '= 0'
[2/4225] Generating qgscontexthelp_texts.cpp
Files "C:/Users/renan.b/QGIS/ms-windows/osgeo4w/build-qgis-dev-x86/src/core/qgscontexthelp_texts.cpp.temp" to "C:/Users/renan.b/QGIS/ms-windows/osgeo4
w/build-qgis-dev-x86/src/core/qgscontexthelp_texts.cpp" are different.
[3/4225] Generating qgsexpression_texts.cpp
Files "C:/Users/renan.b/QGIS/ms-windows/osgeo4w/build-qgis-dev-x86/src/core/qgsexpression_texts.cpp.temp" to "C:/Users/renan.b/QGIS/ms-windows/osgeo4w
/build-qgis-dev-x86/src/core/qgsexpression_texts.cpp" are different.
[6/4225] Generating flex_qgsexpressionlexer.cpp
ninja: build stopped: subcommand failed.
[5/4220] Building CXX object src\native\CMakeFiles\qgis_native.dir\qgsnative.cpp.obj
FAILED: src/native/CMakeFiles/qgis_native.dir/qgsnative.cpp.obj
C:\PROGRA~2\MICROS~1.0\VC\bin\cl.exe   /nologo /TP -DENABLE_TESTS -DNOMINMAX -DQGISDEBUG=1 -DQT_DISABLE_DEPRECATED_BEFORE=0 -DQT_NO_CAST_TO_ASCII -DWI
TH_BINDINGS -DWITH_QTWEBKIT -D_CRT_NONSTDC_NO_WARNINGS -D_CRT_SECURE_NO_WARNINGS -D_USE_MATH_DEFINES -Dqgis_native_EXPORTS -I. -I..\..\..\src\native -
Isrc\native /DWIN32 /D_WINDOWS /W3 /GR /EHsc -DSPATIALITE_VERSION_GE_4_0_0 -DSPATIALITE_VERSION_G_4_1_1 -DSPATIALITE_HAS_INIT_EX  /wd4100  /wd4127  /w
d4190  /wd4231  /wd4244  /wd4251  /wd4275  /wd4505  /wd4510  /wd4512  /wd4610  /wd4706  /wd4800  /MD /ZI /MP /Od /D NDEBUG /D QGISDEBUG -DQGISDEBUG=1
/showIncludes /Fosrc\native\CMakeFiles\qgis_native.dir\qgsnative.cpp.obj /Fdsrc\native\CMakeFiles\qgis_native.dir\ -c C:\Users\renan.b\QGIS\src\native
\qgsnative.cpp
c:\users\renan.b\qgis\src\native\qgsnative.h(34) : error C2065: 'default' : undeclared identifier
c:\users\renan.b\qgis\src\native\qgsnative.h(34) : error C2258: illegal pure syntax, must be '= 0'
[10/4220] Generating moc_qgsdataitem.cpp
ninja: build stopped: subcommand failed.
build failed twice
BUILD ERROR 1: 29/06/2017 14:36:20,63
FINISHED: 29/06/2017 14:36:20,63

Anyone can help me?


texts: ArcGIS Cost Path Analysis. How to extract the area covered by the many possible shortest paths?

texts:
As you can see in the picture ArcGIS Cost Path Analysis produces different results based on chosen direction (blue path direction: start->stop, red path direction: stop->start).

As Cost Path Analysis only works in 8 directions, there are many possible shortest paths in the example above. Some of the shortest paths are located in between the blue and the red path. This got me thinking:

Do you see a possibility to extract the area covered by the many possible shortest paths?

There are large areas of identical costs in this example. Light colors represent low costs, dark colors reperesent high ones.


texts: Different timeout values in print composer and regular map view in QGIS?

texts:
I'm experiencing some issues in the print composer when exporting maps including WMS-layers.
One of the WMS-layers takes a lot of time to load, and I've increased the timeout for network requests in the Network tab (QGIS main window: Settings->Options->Network) to get the WMS-layer to load.
This works when viewing the layers/maps in the print composer as well. If I open the settings window in print composer (Settings->Layout options->Network) I see the same timeout value I entered in the main QGIS window. This makes sense as they most likely point to the same settings.
However, when I'm exporting the map to a PDF from the print composer, the WMS-layer that takes a lot of time to load is not included in the export. I have some vector layers on top of it, and they're included in the export, but not the WMS-layer itself. The other WMS-layers that's in the map view in the print composer are included as normal.
My theory is that the print composer uses a different timeout value when exporting to PDF. Is this true? If yes, where can I increase this value to ensure that the WMS-layer loads properly?
Edit:
To clarify, the problem is that the WMS layer loads properly in QGIS and in the print composer just fine. The issue I'm having is that the WMS layer does not load when I'm exporting the map. The WMS-layer is not included in the exported PDF. (but loads just fine in the print composer)


texts: Draw polygons from self-intersecting lines: Looking for a workaround for inner boundaries / holes

texts:
I use OGR & Python to develop a geodesic polygon buffer. The buffered polygon is generated from a polyline ring. The polyline ring is generated from point coordinates. As the polyline ring is often self-intersecting, this leads to holes in the buffered polygon geometry. For the outer boundary of a self-intersecting line / polygon, applying a Buffer of 0 to the buffered polygon geometry is a widely used workaround to fill holes from self-intersections and to generate a correctly displayed polygon. However, I have difficulties buffering polygons with holes.

A (buffered) polygon with a hole consists of two lines: An outer line from which I want to determine the most outer boundary (using a Buffer of 0) and an inner line, from which i want to determine the most inner boundary. The workaround using a Buffer of 0 only determines the most outer boundary of a self-intersecting polygon. I am thereby looking for a workaround in order to determine the most inner boundary of a self-intersecting polygon (for applying the buffer to the hole).

The polygon I want to buffer:

I applied several buffer widths. The polyline rings and buffer vertices for the respective buffer polygons are displayed here. In order to create a buffered polygon with a hole inside, I need to determine the most outer boundary of the outer line and the most inner boundary of the inner line:

The polylines are self-intersecting and lead to Polygons with small holes (hole in the top right corner of the polygon):

As for the outer boundary, I use a Buffer of 0 in order to determine the most outer boundary of the outer Polyline. At this point, the hole in the top right corner is closed:

However, this approach can not be used for the inner ring of the polygon. I need to determine the most inner boundary of the self-intersecting inner polyline (Lines inside the hole in Image 2). From both, outer line with outer boundary (which I already have) and inner line with inner boundary I want to generate a buffered polygon with a hole inside.

Is there a workaround for this?


texts: Set style from features of layer to hidden?

texts:
Is it possible to set the style to hidden of certain parts of a layer of the map? The API documentation does not seems to have any property for it?

I tried to add the following:

var locatusObject = new ol.layer.Vector ({
  source: locatusSource,
  maxResolution: 14,
  style: styleFunction
});
function styleFunction(feature,resolution) {
style = [new ol.style.Style({
                stroke: changeStrokeStyle(),
                fill: changeFillStyle(),
                //opacity:0,
                //visible:none
                //visibility:hidden
            })];

return style;
}

None seem to work?

EDIT

Oke, it seems to work when I set the variable "style=null"!
But when I click a certain category, my map does not get reloaded?

I use the following code to click and reload my map:

$('#legend ul > li > ul > li').click(function(){
      if($(this).attr('data-hidden')==='false'){console.log('in conditie'); $(this).attr('data-hidden', 'true');$(this).addClass('enabled');}
      else {$(this).attr('data-hidden', 'false');$(this).removeClass('enabled');}

      $('li').each(function(){
        hiddenArray[$(this).attr('data-name')] = $(this).attr('data-hidden');

      })
      map.renderSync();
    });

The map.renderSync() doesnt seem to do anything to my map?


texts: How can I convert the OGC WKT for Icelandic Lambert 1900 to Proj4

texts:
how can I convert OGC WKT to Proj4 string

http://spatialreference.org/ref/epsg/3052/

PROJCS["Reykjavik 1900 / Lambert 1900",GEOGCS["Reykjavik 1900",DATUM["Reykjavik_1900",SPHEROID["Danish 1876",6377019.27,300,AUTHORITY["EPSG","7051"]],TOWGS84[-28,199,5,0,0,0,0],AUTHORITY["EPSG","6657"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4657"],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]


texts: Using QgsDxfExport in Python script leads to crash

texts:
I am trying to export a QgsDxfLayer from QgsVectorLayer, but it crashes at dxf.addLayers().Might I forget something essential ?
...
dxfLayers = []
dxfLayers.append(QgsDxfExport.DxfLayer(vecLayer, -1))

dxf = QgsDxfExport()

# Crashes here :-(
dxf.addLayers(dxfLayers)

I am using QGIS 3.8.
Here is a little bit more code
for layer in QgsProject.instance().mapLayers().values():

    if layer.name() not in LAYERS_OF_INTEREST:
        continue

    info("------- '" + layer.name() + "' ------------------------")

    dxfLayers = []

    # Iterate through all layer features, create vector layer
    # new (vector-)layername, make a vector layer
    # instance and add it to layerList ..
    for feature in layer.getFeatures():
        try:
            # Create our own VectorLayer from feature,
            # make a DXFLayer from it and add it to dxfLayers[]
            vecLayer = MyVectorLayer(feature)
            dxfLayer = QgsDxfExport.DxfLayer(vecLayer.layer, -1)
            dxfLayers.append(dxfLayer)

            if STORE_LAYERS_AS_SHAPEFILES:
                vecLayer.write2File(STORAGE_PATH)

        except Exception as e:
            notHandled += 1
            #print(str(e))
            pass

    if len(dxfLayers) > 0:
        # We have some dxf-layers in our list, so lets create
        # a QgsDxfExport instance and write the dxf-file to disk

        info(". Exporting dxf ...")

        dxf = QgsDxfExport()
        # I don't know why, but it
        # crashes here (FIXME)
        dxf.addLayers(dxfLayers)
        info(". ok 1")

        dxf.setDestinationCrs(layer.crs())
        info(". ok 2")

        # Export Dxf file
        dxfPath = STORAGE_PATH + "/" + layer.name() + ".dxf"
        info("  path '" + dxfPath + "'")
        dxfFile = QFile(dxfPath)
        dxfFile.open(QFile.ReadWrite)
        info(". ok 3")
        dxf.writeToFile(dxfFile,"UTF-8")
        dxfFile.close()
        info(". ok 4")
        exported += 1

info(". DONE, exported=" + str(exported) + " unknown=" + str(notHandled))


texts: "Export to vector file failed. Error: Feature write errors" in QGIS 2.18 but works in 2.4

texts:
I'm having a problem with QGIS 2.18.7 'Las Palmas' which did not occur with version 2.4 (the standard with my OS - Debian Linux).

I have a *.kml file containing the boundaries of UK parliamentary constituencies. It loads into QGIS as a vector layer fine, but if I then try to save  as *.shp, I get the following:

===================

  Export to vector file failed. Error: Feature write errors: Feature
  geometry not imported (OGR error: Field URI of width 255 truncated to
  254.) Feature geometry not imported (OGR error: Field URI of width 255 truncated to 254.) Feature geometry not imported (OGR error: Field URI
  of width 255 truncated to 254.) Feature geometry not imported (OGR
  error: Field URI of width 255 truncated to 254.) Feature geometry not
  imported (OGR error: Field URI of width 255 truncated to 254.) Only
  645 of 650 features written.

================

This doesn't happen with QGIS 2.4.

Why is this happening?


texts: filter by last digit works in points but not in lines

texts:
I have three-point layers T, A1 and A2. T is original points and A1 and A2 are generated by some algorithms. I am trying to compare the result of A1 and A2 with T.

On T, A1 and A2 I can apply a filter at one from a console with this code.

layers = iface.mapCanvas().layers()
for layer in layers:
    layer.setSubsetString('"field_3" % 100 = 1')

After that, I create hub lines T-A1 and T-A2. I want to filter these lines too same as the points, but the filters don't work. I also tried with strings

right("field_3",1) LIKE '1'

in the layer filter but it shows all the lines without filtering.

NOTE: this expression right("field_3",1) works in the label expression


texts: Create a PLSS quarter-quarter grid using QGIS

texts:
Is there a way to create quarter-quarter sections (of the PLSS township/range/section grid) using QGIS? I've already found http://ianbroad.com/creating-quarter-quarter-section-grid-python/ but using ArcGIS is not an option.

Having a way to label all the sections would be great since I need to complete a whole state.

Other resources I've looked at:

Creating grid constrained within polygon in QGIS?
Creating vector grid using QGIS?
QGIS 2: How to create a regular vector grid at fine scale and for large area?
How to create a fishnet grid Shapefile in QGIS?

I am specifically working with data for Montana and cannot find quarter-quarter grids (aka 16 per section).

Also, making a grid is not perfect, I am just checking to see if anyone has attempted to create a tool for QGIS like the one Ian Broad created for ArcGIS. This will not create perfect sections, but gets the job done.

I did find these, but they are only partial files:

http://svc.mt.gov/msl/mtcadastral/
http://www.geocommunicator.gov/GeoComm/lsis_home/home/index.htm


texts: How does raster package convert text to number values on export to tif?

texts:
Background:

spatialpixelsdataframe with @data$column classifying the pixels into n categories
the category column contains text data e.g. "Area A", "Area B" etc and is of factor type
writing that information to GeoTIFF requires that the text names be converted to number values

writeRaster (or writeGDAL? or the GTiff driver?) seems to be able to do this without instruction. So it must decide how to assign numbers to those categories. I'm not sure how that's done. By factor level? Order of encounter as values are written out? I'd like to know for sure, because it affects how one writes lookup tables for the output tifs (Value 2 = "Area A" etc).

I'm aware that one could convert the categories to numbers in a new column before writing that to file, but I still want to know what happens if one doesn't.


texts: Identify feature which is set for editing

texts: While creating custom tool for ArcMap I found problem which couldn't be solved even with try-catch block. When start editing I must to choose which layer I want to edit (if they are from different sources), but if my tool tries to edit other than I've chosen - ArcMap shuts down. How can I identify which layer(s) I can edit?


texts: Finding duplicate shapely polygons (same polygons with different points)

texts:
I am interested in finding out duplicate shapely polygons having different points but covering the same area. For example, let assume I have 2 polygons:
poly1 = POLYGON((4 1, 30 1, 45 1, 45 6, 12 6, 4 1))

poly2 = POLYGON((4 1, 45 1, 45 6, 12 6, 4 1))

These 2 polygons are representing the same area on the plane but are considered different when I do some opertions on them, like these 2 polygons (individually) intersecting with a 3rd polygon.


texts: Finding the UTM coordinates for the start point and the end point of features in a line feature class with ArcPy, 2.7.8, ArcMap 10.8

texts:
I have been trying to create a script to find the start and end UTM coordinates for a line feature class in ArcMap 10.8. The server we are using has Python 2 so my script is in Python 2. I have been trying to use the ConvertCoordinateNotation_management tool with limited success.
Steps:

I first found the coordinates of the start and end of the line. I have used two different ways to do
this. (see code below). adds x and y coordinates for the start and end points to my FC.
I then use the ConvertCoordinateNotation tool to find the UTM coordinates for each feature in the feature class created in step 1. This is done twice to try and get the UTM for the end points and the start points

Unfortunately it is only giving me the coordinates for one of the features in my feature class in the PofC file created and nothing in the PofT feature class created and the UTM coordinate that is added is incorrect. How do I get it to return UTM coordinates for all the  features and for both the start points and end points?
My geographic coordinate system is GCS_North_American_1983 and my projected coordinate system is NAD_1983_BC_Environment_Albers. I know that the location of the one point it is returning is likely because of this but am unsure how to correct it.
# imports
import arcpy

#gdb = raw_input("please input the path of the GDB you are working with. >> ")
gdb = r"W:\Projects\AMD5.gdb" #for testing purposes
arcpy.env.workspace = gdb
#in_fc = raw_input("Please input the fc that requires utm coordinates. >> ").upper().strip()
in_fc = "AMD5_FC" #for testing purposes
#out_fc = raw_input("Please enter the name of the new feature class. >> ").upper().strip()
out_fc = 'test_vertices' #for testing purposes

# set parameter values
# PofC = 'PofC'
# PofT = 'PofT'
# PofC_Y = 'START_Y'
# PofC_X = 'START_X'
# PofT_X = 'END_X'
# PofT_Y = 'END_Y'
# input_format = 'DD_2'
# output_format = 'UTM'

# arcpy.ConvertCoordinateNotation_management(in_fc, PofC, PofC_X, PofC_Y, input_format, output_format)
# arcpy.ConvertCoordinateNotation_management(in_fc, PofT, PofT_X, PofT_Y, input_format, output_format)
# print("The conversion is done.")

#####################################################################################################

arcpy.FeatureVerticesToPoints_management(in_fc, out_fc, "BOTH_ENDS")
PofC = 'PofC'
PofT = 'PofT'
PofC_Y = 'START_Y'
PofC_X = 'START_X'
PofT_X = 'END_X'
PofT_Y = 'END_Y'
input_format = 'DD_2'
output_format = 'UTM'

arcpy.ConvertCoordinateNotation_management(in_fc, PofC, PofC_X, PofC_Y, input_format, output_format)
arcpy.ConvertCoordinateNotation_management(in_fc, PofT, PofT_X, PofT_Y, input_format, output_format)
print ("script complete.")


texts: Selecting Workspace and Datasets from same dialog using ArcObjects?

texts:
I want to select the workspace and the dataset from the same dialog and this is my code for the workspace:

 private IWorkspace BrowseForWorkspace(int hwnd) {

            IGxDialog gxDlg = new GxDialogClass();
            gxDlg.ObjectFilter = new GxFilterWorkspacesClass();
            gxDlg.AllowMultiSelect = false;

            IEnumGxObject workspaces;
            IWorkspace ws = null;

            if (gxDlg.DoModalOpen(hwnd, out workspaces) && workspaces != null) {
                IGxObject gxObj = workspaces.Next();
                if (gxObj is IGxDatabase) {
                    var gxDb = gxObj as IGxDatabase;
                    ws = gxDb.Workspace;
                } else if (gxObj is IGxFolder) {

                    IWorkspaceFactory wsf = new ShapefileWorkspaceFactoryClass();
                    ws = wsf.OpenFromFile(((IGxFile)gxObj).Path, hwnd);

                }
                Marshal.FinalReleaseComObject(workspaces);
            }
            gxDlg.InternalCatalog.Close();
            Marshal.FinalReleaseComObject(gxDlg);
            return ws;

        }

and this is my code for the dataset:

private IDatasetName Browse(int hwnd) {

        IGxDialog gxDlg = new GxDialogClass();
        gxDlg.ObjectFilter = new GxFilterDatasetsAndLayersClass();
        gxDlg.AllowMultiSelect = false;
        IEnumGxObject workspaces;
        IDatasetName datasetName = null;

        if (gxDlg.DoModalOpen(hwnd, out workspaces) && workspaces != null) {

            var gxObj = (IGxDataset)workspaces.Next();
            if (gxObj != null) {
                var dataset = gxObj;
                datasetName = dataset.DatasetName;
            }
            Marshal.FinalReleaseComObject(workspaces);
        }
        gxDlg.InternalCatalog.Close();
        Marshal.FinalReleaseComObject(gxDlg);
        return datasetName;

    }

How do I link them together?


texts: Upgrading Heron to 1.0.5 to use with OpenLayers 2.13.1?

texts:
Can we use Openlayers 2.13.1 instead of openlayers 2.12 with heron 1.0.5 ?

If yes, then what changes we have to make and if no then what are the possibilities of its up gradation in future.

The major problem I have found with openlayers 2.12 is mouse wheel zoom.

It is not as smooth as in openlayers 2.13.1.


texts: How to (re)install gdal in python3 (mac)?

texts:
I am stuck trying to install gdal, running python3.6 on a Mac.

When I run pip install gdal I get errors of

Failed building wheel for gdal

error: command 'gcc' failed with exit status 1

Command "/Users/scottieb/anaconda/bin/python -u -c "import setuptools, tokenize;__file__='/private/var/folders/1v/jks7f3sj1v772x048q
qttmsw0000gn/T/pip-build-8hl9i7bl/gdal/setup.py';f=getattr(tokenize, 'open', open)(__file__);code=f.read().replace('\r\n', '\n');f.c
lose();exec(compile(code, __file__, 'exec'))" install --record /var/folders/1v/jks7f3sj1v772x048qqttmsw0000gn/T/pip-yvs4yb1x-record/
install-record.txt --single-version-externally-managed --compile" failed with error code 1 in /private/var/folders/1v/jks7f3sj1v772x
048qqttmsw0000gn/T/pip-build-8hl9i7bl/gdal/

Some background:

I previously had a working version on mac, using this answer.  But somehow, in trying to install osmnx, I screwed up my libraries and decided to uninstall anaconda and start afresh.  On the latest install I am resetting up and can no longer use networkx to read shapefiles

import networkx as nx
g = nx.read_shp('Road files/geo_export_4d537b7d-a470-4eb9-b147-1d0ea89e6b60.shp')

I tried to install gdal as before (conda install gdal) but get the error:

UnsatisfiableError: The following specifications were found to be in conflict:
  - gdal -> numpy 1.8* -> python 2.6* -> openssl 1.0.1*
  - python 3.6*

Which I suspect means gdal is not currently supported for 3.6?

This is over my head currently, so I'm interested to learn where I went wrong and what exactly is going on to prevent me from retracing my steps here.

Cleaned up my path after advice from gene.
echo $PATH is now just

/Users/scottieb/anaconda/bin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Applications/Postgres.app/Contents/Versions/latest/bin

Even after running conda install gdal, I can not read the shapefile.  I'm getting the error

ImportError: dlopen(/Users/scottieb/anaconda/lib/python3.5/site-packages/osgeo/_gdal.so, 2): Library not loaded: @loader_path/./libgssapi_krb5.2.2.dylib
  Referenced from: /Users/scottieb/anaconda/lib/libgdal.20.dylib
  Reason: image not found

This is pretty common, and the normal solution is to install gdal.  But I've already done that.  Even though the install runs fine, I also can't import gdal while in python, a similar error about libgdal.20.dylib.  It looks like I have 2.0.0 for both gdal and libgdal, but it won't update (ie to 2.1.3, which it looks like gene has below).


texts: Can geoserver limit layer access to clients by credentials?

texts:
I have 2 clients which don't want to host geoserver on their systems and use my server.
Since I can encode their js , is there a way to limit the workspace for each client do other can't call it using wms ?


texts: How to do a one-to-many Spatial Join/Query in QGIS

texts:
Does anyone know if it's possible to do a 'one to many' spatial join based on distance in QGIS?  The "Join Attributes By Location" under ->Vector->Data Management Tools is very limited...

Just to be clear, below is how it's done in ArcGIS.


texts: Linear referencing - Make Route event layer: wrong spacing of points along route?

texts:
I want to simulate the ArcGis tool Generate points along line with linear referencing tools (we don't have the licence for the former).

Our methodology is to

create/use a line
create an excel table with m-values (or the distance from the origin of the line) (see example below)
For our problem, we want equally spaced distances between the points
create a route from the line
Use the Make Route event layer tool to match the m-values to the route in order to create our "points along the line"

The problem is, that any time there is a bend in the line, the spacing of points is wrong (the gap is too large)

See 2 examples below, the right line is our actual business case, the left one I created myself to do a sanity check (digitized the line myself). On the left example you can clearly see that the problem occurs at every bend.

Below you can see the excel file we use to specify the point spacing:

This is the attribute table of the line on the left:

Tool parameters we used:

And the parameters used in the Create Route tool


texts:
I am trying to install the tiger geocoder for postGIS. I am using the 2.1 version of postGIS on Windows 7. I am using this set of instruction to do so.

I am confused by what they mean be edit the paths in the declare_sect column to those that fit your pg, unzip,shp2pgsql, psql, etc path locations.

Here's the code that I am using to do this.

INSERT INTO tiger.loader_platform(os, declare_sect, pgbin, wget, unzip_command, psql, path_sep,
           loader, environ_set_command, county_process_command)
SELECT 'ravimehta', declare_sect, pgbin, wget, unzip_command, psql, path_sep,
       loader, environ_set_command, county_process_command
  FROM tiger.loader_platform
  WHERE os = 'windows';

@RyanDalton - I tried to do what you suggested. However, I got the following error.

  ERROR:  unrecognized configuration parameter "unziptool." I'm not sure
  what I am doing wrong.

Fixed problem by using:

UPDATE tiger.loader_platform SET unzip_command = 'C:\Program Files\7-Zip\7z.exe' WHERE os = 'ravim';


texts: Where are the QGIS Style Manager settings stored?

texts:
I would like to export all the icons defined in QGIS Style Manager to another QGIS installation in another PC.

I've been looking in Windows Registry inside "HKEY_CURRENT_USER\Software\QGIS" but I did not find anything.


texts: Converting Raster Values from Long Integer and Double Values to Short Integer

texts: In order to make my categorical raster data compatible with another model, I need to convert the values from Long Integer and Double to Short Integer. Is there a way to export or convert the files which would do this? I am using ArcMap.


texts: Removing duplicate points in GeoPandas using distance and time


texts: How were walkways found and added to Google Maps or OSM?

texts:
Does Google Maps or OSM use any deep learning or machine learning model to extract walkways from satellite images? Examples of walkways are the green lines in the picture here. I wonder how they were added to the map in scale?
Most papers about pedestrian network extraction I found online are at least 7 years old. I wonder if there are any more recent work in walkway extraction?


texts:
I'm writing a python code to convert X,Y points in excel to shapefile. In this process i have toread xy points from shhet1. What is the step i have to include in my process to read data from sheet1 of an excel work book(97-2003).
I wrote a code as follows..

import arcpy
from arcpy import env
import Xlrd
env.workspace="E:\"
input_table="123.xls\Sheet1"

But its not working.


texts: I want to analyze using world traffic service. There are several points on the map. From these points, the process of determining the area that can go for 15 minutes by car will be made. How can I do this? I use ArcGIS World Traffic Service.


texts: Creating Thiessen polygons gives ArcGIS ERROR 000537 error in converting points to Voronoi Polygons?

texts:
I have a problem with creating Thiessen polygons for a set of points. I am creating them with Create Thiessen Polygons tool in ArcGIS 10.1). I am able to create them in a few steps (for parts of the data) but I would like to create it at once.
The error I get is:

  ERROR 000537 error in converting points to Voronoi Polygons

It doesn't say much, when click at the error in arcgis help I found this message:

  Description
  A TIN surface is created from input points and is used to generate the output Thiessen (Voronoi) polygons.

  Solution
  Check the input point features and make sure that features exist in this particular feature class. As well, make sure that if the extent environment is set, the input point features do not fall outside their coordinate range.

But if it works "partially", where the problem may be?


texts:
Since our ArcGis applications are developed as desktop and as server applications, I'm looking for a way create our SQL-queries in a safe way. The state of the art to mitigate SQL-injection is, to use parametrized queries / prepared statements.

Prepared statements can only be implemented with support of the underlying database and I'm not aware that Esri offers support of those in their interfaces. Parametrized queries could partly be implemented in an external library.

➽ Is anybody aware of a library, which offers parametrized queries for any of the ArcGis products, or does Esri offer support for parametrized queries in one of their products? How did you solve this problem?


texts: Calculate statistics of raster based on plygon borders in google earth engine

texts:
I have a script in GEE that takes shepfile (sometimes with two polygons), classify it to 3 different classes and create one new layer with those 3 classes.
I want to calculte statistic for each polygon with the number of pixels that beling to each class.
I have  a code for this:

 var stats = function(ndviClassImg){

  var num_pos=ndviClassImg.gt(2).selfMask().reduceRegion({
      reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 20
  }).values().get(0);

    var num_neg=ndviClassImg.lt(2).selfMask().reduceRegion({
      reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 20
  }).values().get(0);

    var num_mean=ndviClassImg.eq(2).selfMask().reduceRegion({
      reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 20
  }).values().get(0);

  return ndviClassImg.set('num_pos', num_pos).set('num_neg', num_neg).set('num_mean', num_mean);

};

but when I print the image or the geometry I don't get this data about how many pixels I have in each class.

My end goal: to get the data for each polygon about how many pixels I have in each class (I want also to do more calculation with this later on)


texts: How to join specific attributes of two polyline shapefiles that overlapped but information is stored differently (QGIS)?

texts:
I have just started using QGIS (2.12.3) and I need your help.

Overall, I have two polyline shapefiles with the same geometry (the one layer overlaps the other), but the information is stored in a different way. Furthermore, each layer has different attributes. My purpose is to copy specific data from the one layer to other one.

In particular, let's say that I have a layer called F and a layer called M. The layer F is consisted of sets of many lines, whereas the layer M is consisted of sets of a few lines. See what happens when I select the line1 on layer M and what happens when I select the same line on layer F (Image 1).

This occurs because, the corresponding line1 of layer M in layer F, is part of many polylines (Image 2). Therefore, for one entity in layer M, I have many entities in layer F.

The lines in layer F have a field (integer) and my target is to copy each time (e.g. overlapped lines for line1) only the row that contains the minimum value of this field, into the corresponding line in layer M.

For example, if the field of interest in layer F for the  line 1 has values 1987 (a), 1802 (b), 167 (c), 2989(d), I would like to join to the line1 of layer M the row of layer F that contains the number 167.

I have tried the spatial join, but this tool takes "the attributes of the first located feature". If there was any possibility to take "the attributes of the min value of a specific field" then it would be the solution, I think. Also, I cannot remove the redundant values because with this way I will lose the correct information for other lines.

Any ideas??


texts: ArcGIS API Javascript - FeatureLayer from json: problem with queries and popup

texts:
Hello everybody
I've tried this code proposed by Derek Swingley in JSFIDDLE:
JSFIDDLE-FeatureLayer from JSON
Everything "seems" to work fine, but not always. If you go down to the closer scale (level 21 I suppose) and click several times on different points of the line, sometimes you get a popup with the text "(1 of 3)" in the black section of the popup and sometimes you get the text "Buck Creek 3". Once you press in the arrow keys of the popup, this title doesn't return, you will always have the "(x of 3)" text and, besides, the order is a bit strange. Pressing the right arrow, it passes from "Buck Creek 3" to "Buck Creek 2" and after that to "Buck Creek 3" again. It's needed to press several times the left arrow to get to "Buck Creek 1".
I have built a similar example, replacing lines with points in the json. In this case, I never get the "(1 of 3)" text when the popup has just opened. Now, similarly, if I press the arrow keys of the popup, I always get the "(x of 3)" text. And with points I have same problem with the order.
I think the problem is related to the matter that the feature layer is coming from a JSON, because if I create it from an ArcGIS Server URL I'm not able to reproduce the problem. But the problem is that I need to take the data from the JSON.
Do you know if it's possible to modify this behavior? I prefer having the "(1 of 3)" text from the very beginning and see the elements in the right order.
Thanks in advance,
Santiago


texts: Preventing ModelBuilder model from working on/after particular date?

texts:
I am working in ModelBuilder in ArcGIS Desktop.
I want to stop models from working when the current Time is greater than the project over time. Below snap showing where I am making a mistake please help.


texts:
I am using the ArPy example from esri sample. Script was working fine but then I updated my ArcGIS Server to 10.2.2 and I am getting this error: Token Expired, Code: 498. I have proper credentials as I am able to login to ArcGIS via Manager.

http://resources.arcgis.com/en/help/main/10.1/index.html#/Example_Publish_service_definitions_listed_in_a_text_file/01540000061r000000/


texts: When I press a layer name and select "hide layer" in QField nothing happens, the layer won't hide. I've had experienced it with some layers but not with others. Then I am lost.


texts: GEE-Select pixels whose values equal to a given value in an image

texts:
There is a classified image in which each pixel has a class number.
I want to select those pixels with the same class number.
Then do some statistics such as mean for those selected pixels.
How can I achieve this goal?
Looking forward for your instruction!


texts: Splitting multipart polygons to single part in R?

texts: I am looking for the equivalent of "Convert multipart features to multiple singlepart features. Creates simple polygons and lines" in QGIS and similar in ArcGIS for R. Various analyses create multipart features (in class sf, they are in the class MULTIPOLYGON or MULTILINESTRING) when the desired output is a simple polygons or simple lines (in class sf, they are in the class POLYGON or LINESTRING)


texts: Can you export location details from an image submitted in forms and directly import this information into the forms produced CSV?

texts:
I am hoping to produce a form that can create lines of data including geospatial data in order to plot this information on a map. My form requests an image followed by questions to supply data with the image and this nicely populates into a CSV.
What I would now like to do is grab location data from the image and input into my CSV so that I can plot the information at the location in which the image was taken?
If there is a better method than using forms for this task I am happy with that, however I need the data input to be user friendly.


texts: Output added to display & symbology applied only when output is a Model Parameter

texts:
When I run my model the output is not added to display even though the Add to Display option is checked. Also, symbology applied from a layer I specified is not showing. Checking the Model Parameter option solves both of these issues. However, since I am going to publish the model as a geoprocessing service I don't want the user to have the option to mess with the name and location of the output.

It is not the first time this is happening nor am I the only one it happens to so I don't think there is a problem in my model.

I am using ArcMap 10.4.1.


texts: How to workaround large mosaic process which is failing?

texts:
I need to mosaic about 550Gb of tif imagery together and the software I have tried keeps failing. The area has been split into zones so that the smallest has approx 200 tiles.

I have used the latest versions of ERDAS (Imagine and Mapper), ArcINFO and Global Mapper on a 3.30 gigahertz Intel Xeon E31245, DELL, 16GB RAM, 64-bit Win 7 Professional. Mullti-core (4 total), Hyper-threaded (8 total) machine. My C has 700GB free and D has 1.5TB.

I am looking into using Grass (never have before) but i.image.mosaic only seems to only handle 4 files...some of mine have 600 tiles. Any other options or opensource software to try?

Sorry should add that we can't use a mosaic dataset (or equivalent in other software) as we need to create zones with defined no-data areas as ecw's so that they can be opened in any GIS software and combined with lower resolution/older data when new data doesn't exist seamlessly.

 An example of how some mosaiced files look in different sofware. Global Mapper/ERDAS are fine but it's not correct in arcgis.

--- OLDER INFO---

Sorry for the rough drawing. So having the colored areas as 5 zones will minimize the no data areas in the larger AOI.

In arcgis the code is as follows (this is run as a model and not in python as I can't get it to take the tifList input).

arcpy.MosaicToNewRaster_management(tifList+";" +mask,RootOutput,"Tile1.tif","PROJCS['GDA_1994_MGA_Zone_55',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',147.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]","16_BIT_UNSIGNED","0.5","3","MAXIMUM","#")

# Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
# The following inputs are layers or table views: "test2"

arcpy.CopyRaster_management(OutputFile,RootOutput+"Tile1b.tif","#","256","256","NONE","NONE","16_BIT_UNSIGNED")

where tifList should be read in from a csv file but this didn't work in python so I am running the above in a model instead...

I have 1.5TB+ free space on my drive but the process crashes with a 9999 error.

Would even 100 tiles process? -i.e should we look at breaking the zones up further?


texts: projection for high accuracy small areas

texts:
I need to convert (lat,long,height)points into a cartesian system for an area as small as 40kmx40km or even 10kmx10km. I'm currently using the UTM projection, but there are some issues regarding accuracy and distortion. The good thing with UTM is that GDAL library easily do the math via python coding.

My professor told me 'plano topográfico local' (in Portuguese) would be a good option. The bad thing is that I have no idea of the EPSG code for this system in order to use GDAL. Any ideas?


texts: Determining which polygons only intersect once with the polygons of another shapefile

texts: I have two shapefiles, one for counties and one for congressional districts. I need to determine which counties are in more than one congressional district. I would think I could use postgis and use some combination of ST_intersects and ST_touches, or ST_overlaps, unless there is an easy option in Arc or QGIS. What makes this more difficult is that I need to do this for 100 years back so the further back you go the less data there is, so for example, I have a county shapefile for 1970 that has a resolution of 1:250,000 and a congressional district shapefile for 1970 that has a resolution of 1:500,000, and I can't find shapefiles that match resolution, so how would I work around that problem as well. Just wondering if anybody has any pointers for either of those problems.


texts: Interview questions for hiring GIS Analyst?

texts:
What are some suitable interview questions to ask candidates for a GIS Analyst position?

I am looking for techniques for interviewing analysts at varied experience levels (we are currently looking at hiring for entry and senior-level positions).

In the past, I've asked questions about their most recent projects and what sites/newsletters they use to keep up with the industry.


texts:
I need to make an esriRequest without filtering the results, like this:

require([
    "esri/config",
    "esri/request"
], function(esriConfig, esriRequest) {
    esriConfig.request.proxyUrl = "/proxy/";
    esriRequest(
        'http://.../MapServer/0/query',
        {
            query: {
                f: "json",
                where: "(1=1)",
                outFields: '*'
            },
            responseType: 'json'
        }
    ).then(function(response) {

    });
});

My problem is that features.geometry data is really big, and it takes several seconds before end-users see the result.
Is there a way to only receive feature.attributes data?


texts: How is a GeoTIFF geographically referenced if it has no sidecar files?

texts:
I am currently downloading MODIS burned area fire products to use for a work project. One of these products is a GeoTIFF which indicates the day (Julian date) of the year when an area experienced fire.

These GeoTIFF files are available for download by themselves with an extension of .tif and without any sidecar files (e.g. .tfw, etc.). When I open the .tif file in ArcPro, it is correctly georeferenced - how is this possible without a sidecar file which contains georeferencing information?

I merely ask this as to me, conceptually, a .tif file is analogous to any other simple image file (e.g. .jpg) and I am curious how it can be natively georeferenced?


texts: How to align QGIS legend items to the right?

texts:
Is there a way to align QGIS map composer legend items to the right, for right-to-left languages?
There is only a title alignment option.
Furthermore, I found this ticket from 4 years ago regarding same issue: trac.osgeo.org/qgis/ticket/3239

My current workaround is to add spaces in item name, but that only aligns the text, not the symbol obviously.


texts:
I am working on a forest project. I have one central point (I know coordinates of it), than I have azimuth and distance of trees from this central point.

How can I find coordinates of trees?

p.s. I am working on QGIS.


texts: Displaying Distance in Miles for Elevation Profile using Profile Tool in QGIS

texts:
I am attempting to generate an elevation profile in QGIS 3.2.2, with elevation in feet, and distance in miles, using Profile Tool.  I have managed to convert my DEM files to feet, and the elevation of the path displays correctly in feet in the Profile Tool window.  The distance continues to display in meters.

I have set my project preferences to use Miles for distance measurement, but no matter what I do, what setting I click, I cannot get the distance on the x axis to display in feet in Profile Tool.

I've also tried the qProf plugin with the same results; elevation displays in feet, distance in meters.

I feel like I'm overlooking something simple.  Is there an easy way to change the distance units, either in the data, in QGIS itself, or using the Profile Tool?


texts: how to bind custom marker controls to map : leaflet

texts:
I am working on map control app using leaflet plugin ....

I added a custom button for marker in side toolbar....

Target

I want to put marker on map when i click on that button....in other words i want to enable marker control when i click on this button

Problem

I am trying to bind the marker events with that button and as usual i face failure in doing that,,,,(please consider i am using the poly area functions but want to keep markers separate)

In script :

    var greenIcon = L.icon({
        iconUrl: 'assets/img/marker-icon.png',
        shadowUrl: 'assets/img/marker-shadow.png',
      });
    map.on('click', putMarker);
    function putMarker(e) {
        L.marker(e.latlng, { icon: greenIcon }).addTo(map);
    }
    //custom button for marker
    L.easyButton('glyphicon-facetime-video', function () {
    }).addTo(map);


texts: Arcobjects in Python: assign interface to objects with feature (shape) geometry

texts:
I am trying to use ArcObject in Python. I do have problem with assigning an interface (i.e. IZ.InterpolateZsBetween) to my object with geometry data, code is:

from Comtypes.client import CreateObject, GetModule
import arcgisscripting

gp = arcgisscripting.create(9.3)
gp.setproduct("ArcInfo")

esriSystem = GetModule("C:/Program Files (x86)/ArcGIS/com/esriSystem.olb")
esriGeometry = GetModule("C:/Program Files (x86)/ArcGIS/com/esriGeometry.olb")
esriDataSourcesGDB = GetModule("C:/Program Files (x86)/ArcGIS/com/esriDataSourcesGDB.olb")
esriGeoDatabase = GetModule("C:/Program Files (x86)/ArcGIS/com/esriGeoDatabase.olb")

sPath ="c:/temp/test.gdb"
pWSF = CreateObject(esriDataSourcesGDB.FileGDBWorkspaceFactory,\
                    interface=esriGeoDatabase.IWorkspaceFactory)
pWS = pWSF.OpenFromFile(sPath, 0)
pFWS = pWS.QueryInterface(esriGeoDatabase.IFeatureWorkspace)
pFClass = pFWS.OpenFeatureClass("test_z")
pFCursor = pFClass.Search(None, True)
pFeat = pFCursor.NextFeature()
pIZ = pFeat.QueryInferface(esriGeometry.IZ2)

Why this is not working?
The output is:

 Traceback (most recent call last):
  File "<pyshell#175>", line 1, in <module>
    pIZ = pFeat.QueryInterface(esriGeometry.IZ2)
   File "C:\Python25\lib\site-packages\comtypes\__init__.py", line 1078, in QueryInterface
 self.__com_QueryInterface(byref(iid), byref(p))
 COMError: (-2147467262, 'Taki interfejs nie jest obs\xb3ugiwany.', (None, None, None, 0,
 None))

So I tried:

pShape = pFeat.ShapeCopy
pIZ = pShape.QueryInterface(esriGeometry.IZ2)
pShape.InterpolateZsBetween(0, 0, 0, 4)

And this is working.
But how to stick it back and save it in my feature?
I am learning Python and ArcObjects, so please be understandful.
Any help will be appreciated.


texts: How to load a TIFF image in Conda Python through QGIS package

texts:
I installed QGIS package in the Python environment, but while I load a TIFF image, it always fails, and no error reported.

from PyQt5.QtCore import QFileInfo
from qgis.core import *
from osgeo import gdal

if __name__ == '__main__':
    QgsApplication.setPrefixPath\
        ("D:/software/installdic2/Anaconda3/envs/py36/Library", True)
    qgs = QgsApplication([], False)
    qgs.initQgis()
    gdal.AllRegister()
    my_project = QgsProject()
    my_project.write("D:/qgistest/my_project.qgs")

    my_pro = QgsProject()
    my_pro.read("D:/qgistest/my_project.qgs")

    filename = "D:/qgistest/000E00N.tiff.tif"
    file_info = QFileInfo(filename)
    base_name = file_info.baseName()

    layer = QgsRasterLayer(filename, base_name)

    if not layer.isValid():
        print("Layer failed to load")
    print(1)


texts: How to add street view in google maps for some specific area?

texts: I need to know how to add a street view of some specific area on google maps or google earth. Is that possible?


texts: Input multiple KML files into Leaflet API

texts:
As we know, since the WIndy Leaflet-kml plugin has been developed we are able to input our KML files into the Leaflet map.

However I have one query regarding this:

How we can input a huge (at least several).kml files into this map?

Comparing it i.e with Google MyMaps, where we can place 8 single .kml layers (like 8 single .kml files). Could we do it the same in Leaflet using this plugin?

I tried this code, which works:

// Load kml file
fetch('lotA.kml')
      .then( res => res.text() )
      .then( kmltext => {

            // Create new kml overlay
            parser = new DOMParser();
            kml = parser.parseFromString(kmltext,"text/xml");

            console.log(kml)

            const track = new L.KML(kml)
            map.addLayer(track)

            // Adjust map to show the kml
            const bounds = track.getBounds()
            map.fitBounds( bounds )

      })

fetch('lotB.kml')
      .then( res => res.text() )
      .then( kmltext => {

            // Create new kml overlay
            parser = new DOMParser();
            kml = parser.parseFromString(kmltext,"text/xml");

            console.log(kml)

            const track = new L.KML(kml)
            map.addLayer(track)

            // Adjust map to show the kml
            const bounds = track.getBounds()
            map.fitBounds( bounds )

      })

... 2 files more, and last one:

fetch('LotD.kml')
      .then( res => res.text() )
      .then( kmltext => {

            // Create new kml overlay
            parser = new DOMParser();
            kml = parser.parseFromString(kmltext,"text/xml");

            console.log(kml)

            const track = new L.KML(kml)
            map.addLayer(track)

            // Adjust map to show the kml
            const bounds = track.getBounds()
            map.fitBounds( bounds )

      })

Despite the code, which is valid, it can look awkward, when inpuit more files.~

I believe, it must be some solution with the fetch() reference. I was searching some info, but everywhere the fetch() option is related to some url address, on top of that, a single one:

https://developers.google.com/web/updates/2015/03/introduction-to-fetch

https://github.github.io/fetch/

https://scotch.io/tutorials/how-to-use-the-javascript-fetch-api-to-get-data

I was trying to do sth like this:

fetch(['lotA.kml'].['lotB.kml']), but this way is not good.

Do you know how to refactorize this code to make it more confined?


texts:
I am interested in the terra R package as it appears to be a successor to the raster package that is being more actively developed than the raster package. But it seems like some features are missing.
For example I have an aggregated National Landcover Dataset raster in .GRD format.
If I call x <- raster::raster('nlcd_agg.grd') the resulting object has a slot x@data@attributes with the raster attribute table that includes the names of the land cover types corresponding to the integer values in the raster.
However if I call x <- terra::rast('nlcd_agg.grd') I cannot find that attribute table in the object nor do I see how to associate the object x with the attribute table of the raster.
I would like to be able to load the raster and have the associated attribute table with land cover type names and default plotting colors load as well. I understand that .GRD is a native format for the raster package but is there some way to achieve this with the terra or stars packages?
Here are two small files that can replicate the behavior (nlcd_agg.grd and the associated nlcd_agg.gri):

nlcd_agg.grd (Dropbox)
nlcd_agg.gri (Dropbox)


texts: Creating Flooding depth Rasters at every feet

texts:
I am explaining the issue I am facing to automate the process.

I have created a DEM from LIDAR and the Asset elevation is incorporated into the surface.

Purpose of this task: Once I get the new asset Elevation, then I need to incorporate those asset elevations in the previously created LiDAR surface.

We need to automate the below mentioned process so that every time the surface is modified, the flood depth raster’s from 1 feet to 12 feet will automatically generate after running the script.

Here you will observe that Raster surface elevation is indirectly proportional to the Flood depth rasters.

One more thing you will observe that the flooding depth color is same for each feet.

Please let me know if you need more explanation.

How to set up the raster classification in this case?

Precipitation/Stormwater color used:

Setting Exclusion properties

Flooding depth rasters at 12 ft and 5ft


texts: How to get latitude and longitude coordinates from district names in India?

texts: I have created a CSV file with the census data for India. One of the columns is District name in the CSV. However, since I do not have lat and long in the CSV file, QGIS is unable to accept and add it as a delimited layer. Is there an easy way to get lat and long coordinates from the district names?


texts: Finish inserting new line-feature (PyQgis + GUI)?

texts:
I'm writing a small graphical plugin where users should be able to add new (multi-)line features to an existing layer.

This is what I've already done:

I have two buttons, one to trigger the insert, the other to stop the inserting when the user is satisfied with his line-feature.
The first button works properly. It selects the associated layer, starts editing it and triggers the Add Feature Button:

self.iface.setActiveLayer(scenario_layer)
scenario_layer.startEditing()
self.iface.actionAddFeature().trigger()

Then the user is able to just click in the map to add vertices.

My problem:

I want the user being able to insert the line using a button, instead of a right click. I cannot figure out how to trigger that in PyQgis.
I still can commitChanges() to the layer, but this won't insert the line to the layer.

Any suggestions?