texts: Find tangent point on circle furthest east or west

texts:
There are few ways of articulating my question. It starts with a known circle on a sphere, i.e. a point-radius, and one that does not cover a pole.  I would like to calculate the point furthest west on this circle and the point furthest east. These points are the tangent points of lines of longitude which encircle the poles and are vertical.  These points are opposite each other along the vertical axis so there is one hard problem to solve from which the other is trivially determined.  This is highly related to calculating the minimum bounding lat-lon box of a circle.  I know how to determine the longitude delta from the center to the edges, which is:

asin( sin( dist_rad / cos(center_lat_rad) ) )

whereas dist_rad is the circle's radial distance is radians, and center_lat_rad is the latitude point of the circle's center in radians.  This equation was elusive to me until I discovered it with a wonderful explanation here.  It contains a wonderful illustration:

.

The illustration shows T1 and T2 which are the points I am trying to calculate. What is left to figure out is at what latitude the eastern and western boundary meridians intersect this circle.  If the circle is at the equator, then it's zero, and if it's north of the equator, then this latitude is going to be higher than the circle's center latitude.


texts: Access point layer representation in QGIS with pyqgis

texts:
I am exporting the style of a categorized point layer with QGIS:

unfortunately my code only extract the size of the symbol in the definition of the properties of the layer but not the representation in the map...
here is a working code:

canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()
i = allLayers[0]
iter = i.getFeatures()
provider = i.dataProvider()
attrvalindex = provider.fieldNameIndex(i.rendererV2().classAttribute())
categories = i.rendererV2().categories()
color_str = []
radius_str = []
transp_str2 = []
transp_str = str(1 - ( float(i.layerTransparency()) / 100 ) )
for feat in iter:
    fid = feat.id()
    attribute_map = feat.attributes()
    catindex = i.rendererV2().categoryIndexForValue(str(attribute_map[attrvalindex]))
    if catindex != -1:
        color_str.append(str(categories[catindex].symbol().color().name()))
        radius_str.append(str(categories[catindex].symbol().size() * 2))
        transp_str2.append(str(categories[catindex].symbol().alpha()))
    else:
        color_str.append('#FF00FF')
        radius_str.append('4')
        transp_str2.append('1')


texts: Exporting layer to shapefile using PyQGIS?

texts:
I´m using python in QGIS 2.2 and struggling with an error trying this:

# create layer from text
_vlayer = QgsVectorLayer(_ur1, "raw", "delimitedtext")

works fine
but:

# export layer to shape
_writer = QgsVectorFileWriter.writeAsVectorFormat\
(_vlayer,"hoppla.shp","utf-8",None,"ESRI Shapefile")

... ends up wiht error code "3" - according to documentation:"ErrCreateLayer".


texts: Auto-generating map images with different styles in QGIS?

texts:
I have got a simple map in QGIS, containing only 1 layer which consists of 10+ features. Now I need to (programmatically) generate about 200 different JPEG map images of the same area. Each map should have the very same map extents and all features visible. However, the features should have different set of colors in each generated map. In detail, I have a CSV file containing the categorized colors for each feature (1 line for each image, 1 column and color per feature), something like this:

"A","B","C","D","E","F" (header line, 6 columns in this file)
"#000000","#FF0000","#00FF00","#0000FF","#0000FF","#808080"
"#000000","#FF0000","#FF0000","#00FF00","#00FF00","#808080"
"#FF0000","#FF0000","#000000","#800080","#FFFFFF","#FFFF00"
...

My current plan is to use Perl or similar similar to create 200 different QML files based on this CSV file. Then I plan to write a Python script to

Loop through QML files
Open each QML file, one at a time, and apply layer style from it
Open map composer and save raster image file
Close map composer and jump into the next QML file

I am now asking for any input concerning:

Can this be done within QGIS + Python?
Any better ways to achieve my goal?
Hints on any existing source code or plugin to achieve this?


texts: QGIS Naming in Legend Based on Multiple Columns

texts:
Is it possible to have the symbol name in QGIS based on two different columns in the attribute table.
I have this right now and I want all the entries to be like the first A. I could manually type them in but I have over 120 entries. Is there a way with an expression to do this for every entry?


texts:
My program needs a numpy array for processing, so I convert GEE image into numpy array using this code:

# Define an image.
img1 = ee.Image( 'LANDSAT/LC08/C01/T1_SR/LC08_038029_20180810') \
  .select(['B2','B3','B4', 'B5', 'B6'])
new=img1.toArray().toArray(1)
# Define an area of interest.
aoi = ee.Geometry.Polygon(
  [[[-110.8, 44.7],
    [-110.8, 44.6],
    [-110.6, 44.6],
    [-110.6, 44.7]]], None, False)

# Get 2-d pixel array for AOI - returns feature with 2-D pixel array as property per band.
band_arrs = img1.sampleRectangle(region=aoi)

# Get individual band arrays.
band_arr_b2 = band_arrs.get('B2')
band_arr_b3 = band_arrs.get('B3')
band_arr_b4 = band_arrs.get('B4')

# Transfer the arrays from server to client and cast as np array.
np_arr_b2 = np.array(band_arr_b2.getInfo())
np_arr_b3 = np.array(band_arr_b3.getInfo())
np_arr_b4 = np.array(band_arr_b4.getInfo())
# Expand the dimensions of the images so they can be concatenated into 3-D.
np_arr_b2 = np.expand_dims(np_arr_b2, 2)
np_arr_b3 = np.expand_dims(np_arr_b3, 2)
np_arr_b4 = np.expand_dims(np_arr_b4, 2)
# Stack the individual bands to make a 3-D array.
rgb_img = np.concatenate((  np_arr_b4,np_arr_b3,np_arr_b2), 2)

Now can I go back from numpy array to GEE image?


texts:
I tried to use QgsProject.instance().addMapLayer(layer) on QGIS 3.4 but I didn't succeed.

Example:

layer = QgsVectorLayer(r'C:\CURSO_PYQGIS\CAPAS\urbano.shp','urbano','ogr ')
QgsProject.instance().addMapLayer(layer)

But the layer didn't show up, and it didn't give me an error on Python console. I don't understand it, someone knows why?


texts: Can't import QGIS library using "C:\OSGeo4W64\bin\python.exe"

texts:
I have ArcGIS and QGIS on my Windows 7 system (installed with OSGeo4W64).
Why I can't import qgis library using "C:\OSGeo4W64\bin\python.exe" ?

How to fix that?
P.S. Importing QGIS from OSGeo4W Shell work well (but it's not what I want).


texts: How to convert building outlines to Polygons?

texts:
I have line vector layer which contains contours of buildings
in the form of separate segments. I want to convert lines to polygons (QGiS 2.17) but it doesn't work.

I'm trying to solve this problem as follows: first off by using algorithm dissolve for all objects, and afterwards by using algorithm multipart to singlepart. Unfortunately, dissolving takes a long time.

How can I do this faster and easier?


texts: Issues loading Layers with Cyrillic Characters to Geoserver 2.12

texts:
I am running Windows 7. With a  Pure installation of Geoserver 2.12.1 (2.12.0 is the same).

Trying to upload local layer with cyrrilic name:

Pure installation of Geoserver 2.10.0 (2.11.4 is the same). The similar uploading:

First image is showing that geoserver after 2.12 version has troubles with cyrrilic layer's names.
Second is showing that if version is under 2.12, problem is gone.

Same data_dir, layers files and pure installations with default settings.

-Dfile.encoding=UTF8 config is useless.


texts:
I am trying to use the 'kernel interpolation with barriers' tool in the Geostatistical Analyst toolbox of ArcGIS 10.3. My data consists of 600 marine animal relocations, and I am trying to calculate the kernel density estimate using a polygon shapefile of the coastline as a barrier (the animals do not use land).

My relocations all have the same weight (i.e., the z- value for all points is 1).

I initially tried the tool by setting the function to 'gaussian' and the 'order of polynomial' to 0. this was based on the parameters used in the paper "Three-dimensional foraging habitat use and niche partitioning in two sympatric seabird species, Phalacrocorax auritus and P. penicillatus" by Adam G. Peck-Richardson et al (2017).

After running for 4.5 hrs, the output was completely useless, where all the cells had a value of 1.

I don't see a lot of documentation for this tool, and I don't understand why I got an output with a uniform value. Any insight?


texts: Converting File Geodatabase Raster to TIFF

texts:
I am trying to export a raster file, saved as File Geodatabase Raster (FGDBR), as a .tif
The usual tools have not worked for me:

data>export feature (resulting error: failed to create raster data set)
raster to other format (to TIFF) (resulting in no file being created)

How do I export this File Geodatabase Raster as a .tif?


texts: How can I find "UUID" in "gdb_items" ? In "gdb_items", "UUID" values are missing

texts:
I need to find "UUID" values in "gdb_items", can I find the values via relationship with other tables?

Note: I have other values in "gdb_items", all fields are fine except the "UUID" field. I just need "UUID" values.


texts: Generating SLD dynamically using PHP in GeoSever

texts:
I am fairly new to GeoSever and I have a question before I launch into something bigger than me: if I replace the XML code in the SLD style sheet inside GeoSever with some PHP which generates the SLD dynamically (i.e. from a database query), will it work?

How does GeoSever handle the SLD stored as style?


texts: Displaying rainbow colors for multiple polygon layers in ArcMap?

texts:
I am working with ArcMap and I have a series (20) layers - polygon lines in a group. I would like to ask you a way to color (not manually) all the lines with a rainbow color palette. The lines are map contour lines (height-mountain), and I would like only to color the lines, not the space between them.

Is any auto way to do that?


texts: Extracting Multivalues to Points from NETCDF4 Raster using ArcPy seems imprecise?

texts:
I'm getting some strange results in a multivalues to points extraction from a NETCDF4 layer from the MERRA reanalysis product. In particular, it seems like the values being assigned are way off in some areas, as you can see from the different colours of the points and the raster, which should be the same since they share the same scale:

It seems like its not a reprojection problem, since the points seem to be where they should be, judging from the contour map of Peru also included.

The code I'm using to extract the multivalues to points (and I'm looping over a long bunch of rasters) is the following:

import os, sys
import Tkinter, tkFileDialog
import glob
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Ask user to select a path to Working directory (i.e. where the raw .nc4 files are)
root = Tkinter.Tk()
root.withdraw()
working_dir =    tkFileDialog.askdirectory(parent=root,initialdir="/",title='Please select a directory')

# Ask user to select CSV file with ENAHO X/Y points
root = Tkinter.Tk()
root.withdraw()
CSV_file = tkFileDialog.askopenfilename(parent=root,initialdir="/",title='Please select CSV file')

# Set the arcgis workspace
arcpy.env.workspace = working_dir
print "Your working directory is:"
print working_dir

####################################################
# Make a XY layer from the CSV
try:
        # Set variables
        x_coords = "longitud"
        y_coords = "latitud"
        out_layer = "enahoXY"
        # 4326 is the code for GCS_WGS_1984
        sr = arcpy.SpatialReference(4326)

        # Make the XY layer
        arcpy.MakeXYEventLayer_management(CSV_file, x_coords, y_coords, out_layer, "", "")

except:
        # If an error occurs print message to screen
        print arcpy.GetMessages()

# Save as .shp
shapefile = os.path.join(working_dir,'/results.shp')
arcpy.CopyFeatures_management(out_layer, shapefile)

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

# Loop around files in directory, but only if they have a .nc extension
files_in_dir = os.listdir(working_dir)
print "Current file processing:"
for f in files_in_dir:
        if f.endswith('.nc4'):

                # Get a part of file that is unique (eg. A2007033)
                parts = f.split('.',2)
                header = parts[2]
                data1 = "tmax"
                extension = ".tif"

                # Program extraction
                temp = "T2MMAX" # maximum daily temperature at 2m (in K)
                x_dimension = "lon"
                y_dimension = "lat"
                band_dimension = "#"
                dimension = "#"
                valueSelectionMethod = "BY_VALUE"

                # Extract temperature
                subset = header + data1 + extension
                print subset
                arcpy.MakeNetCDFRasterLayer_md(f, temp, x_dimension, y_dimension, subset, band_dimension, dimension, valueSelectionMethod)
                # Extract values from Rasters - the column takes the name from the subset
                ExtractMultiValuesToPoints(shapefile, subset, "NONE")

Any ideas what could be the problem?

I'm pretty sure I've entered the correct CRS for the shapefile in the code, but ArcMap tells me the "results.shp" file has no geographical references, which seems strange.

I'm not using the h5py nor the netCDF4 packages for python in the code above.


texts: Set label font size for layer via Python

texts:
I have some ArcMap addin that creates new layer and add it into table of content. I need to set up labels for that layer.

I have made some script for labeling layer:

mxd = arcpy.mapping.MapDocument("CURRENT")
data_frame = mxd.activeDataFrame

layer_divide = arcpy.mapping.ListLayers(mxd, "layer", data_frame)[0]
if layer_divide.supports("LABELCLASSES"):
    for lblclass in layer_divide.labelClasses:
        lblclass.expression = "[OBJECTID]"
        lblclass.showClassLabels = True

layer_divide.showLabels = True
arcpy.RefreshActiveView()

This add labels, but I don't know how to change font size of that labels.


texts: Visualizing custom geo data on WFS featurecollection

texts:
I have some data that I want to visualize over on an interactive map. Basically I want to highlight certain buildings based on the data that I have. My data consists of coordinates and I need to match each coordinate to a building and then highlight that building. Just to illustrate  what I want to achieve, here is a picture I drew in paint.

I am getting the outline of the buildings by using the WFS GetFeature functionality from a publicly available geoserver. I was thinking that I could use the bbox parameter to match each of my coordinates with one of the buildings and then simply draw colored polygons on matching buildings. But considering I have potentially thousands of coordinate pairs this doesn't seem like the right thing to do.

Could someone suggest a good way to achieve what I am looking for? I am quite new to this and have spent quite some time figuring this out, but I can't figure out the best way to do this (One idea I did have is to create my own feature layer using WFS).


texts:
I am using QGIS version 2.18.19. I have created a shapefile using a script and its attribute table is empty as shown in below image.

I am adding an Excel file. It includes a table as the input to the script appeared as follows:

I want to join this table in to my above mention shapefile. I want my shapefile attribute table to be visible as shown below after joining Excel table in to my shapefile.

How can I do this using QGIS?


texts: Lighting Study for Dynamic Sign Measured In Footcandles

texts:
I was assigned to create the attached lighting study to show how bright, and at what distance, an LED-style public sign would illuminate the surrounding area.  The numbers and data came from a chart I was given by the producer of the LED screen.  I had to hand draw it in Google Earth, using angle measurements and distances with lines and circles.  Then I had to manually label each intersection for angle in degrees, and distances with the Measure tool.  The the yellow line is 0° (Y) from the axis of the LED screen (X).  Also had to label the level of luminescence in "footcandles" which are labeled in white on the attached image.

Is there a way I can set up a model in QGIS to automate most of this?

Not only would it save time for when I have to do this again, the results would also be more accurate.  Though I have a lot of experience with ArcMap desktop, modeling was not part of my former position. I have some limited experience with it however.  I also don't have access to ArcMap desktop in my new position.  I currently am using QGIS and Google Earth for mapping tasks.


texts:
Dedicated page on use of leaflet in R provides a very useful example on generating random latitude and longitude values for Central Europe:

cbind(rnorm(40) * 2 + 13, rnorm(40) + 48)

which results in neatly distributed points:

(Source: Leaflet official page)

I'm interested in generating a similar set of randomly generated points by, broadly, corresponding to the UK. For instance, how can I use the latitude and longitude map available through Maps of World to arrive at the X, Y and Z values in the formula below so the generated dots would, roughly, fell across the British Isles?

cbind(rnorm(40) * X + Y, rnorm(40) + Z)

I'm interested in a rough approximation, if some points will fall in the see or outside the UK that's fine.


texts: Why is my Leaflet JS App Canvas moving around?

texts:
I have an Leaflet JS app, and when a user clicks on the canvas, the entire world is able to move causing the unsightly gray background to be shown.

http://jsfiddle.net/nos5gbjh/


texts: Working with varied slope impacts for a cost surface

texts:
I'm working in ArcMap and trying to create a Cost surface for movement through a terrain. I have values (speed in km/hr) for each surface type (tarmac, grass, forest etc).

The problem is moving into the 3rd dimension, z values and slope, created from a DEM dataset. The issue is that the impact which slope has on speed is different moving up slope (every 1% slows by 0.65km/hr) to moving down slope (every 1% increases speed by 0.35km/hr).

Has anyone come across an easy way to deal with this issue?


texts: Marking distances with QChainage Tool and Marker Lines (they differ)

texts:
I have an issue with the markers created by the QChainage Tool and Marker Lines created in Layer Properties (under Symbology). I was trying to use the Marker Lines to create unlabeled intervals between the QChainage markers. The Marker Lines line up fine with the QChainage markers as long as I am zoomed out enough to see the origin. However, when I zoom in and lose sight of the origin, they no longer line up. See the zoomed-out image and zoomed-in image below:

Is there any way to prevent this problem, or am I trying to do something that is just not possible using Marker Lines?

I am using QGIS 3.10.4.


texts:
I know the “combine” tool exists within the spatial analyst toolbox for ArcGIS Pro, but I seem to have run into an upper limit on what it can handle.
Essentially I’ve got a 200 band raster and need to know the number of unique conditions that actually occur (what the unique conditions themselves are is unimportant for the purposes of this question). The upper limit will obviously be the product of the number of classes for each map, but many of these theoretically possible unique conditions won’t occur in reality.
ArcGIS Pro 2.7 crashes about half way through processing using “combine”, presumably because the number of unique conditions is too large for it to handle (I suspect it is potentially trillions...or worse). I have also tried splitting the 200 bands into individual rasters for input to the combine tool and the same problem persists. Combine works on a much smaller number of inputs - i.e. 45 bands creates ~10,000,000 unique conditions, but struggles much beyond that.
Is there an upper limit on the combine tool and, if so, can it be worked around?


texts: Problems with layer labeling QGIS standalone C++ application

texts:
mypLayer->setCustomProperty("labeling", "pal");
    mypLayer->setCustomProperty("labeling/enabled", "True");
    mypLayer->setCustomProperty("labeling/fontFamily", "Cantarell");
    mypLayer->setCustomProperty("labeling/fontSize", "11");
    mypLayer->setCustomProperty("labeling/fieldName", "DRVAL1");
    mypLayer->setCustomProperty("labeling/placement", "1");
    mpMapCanvas->refresh();

Didn't work for me with QGIS 2.14. At the same time this code works fine at Python console in QGIS Desktop application on CentOS 7.5 x64.

layer = iface.activeLayer()
layer.setCustomProperty("labeling", "pal")
layer.setCustomProperty("labeling/enabled", True)
layer.setCustomProperty("labeling/fontFamily", "Cantarell")
layer.setCustomProperty("labeling/fontSize", "11")
layer.setCustomProperty("labeling/fieldName", "DRVAL1")
layer.setCustomProperty("labeling/placement", "1")
iface.mapCanvas().refresh()

Any suggestions?


texts: Opening Notepad text editor from ArcPy/Python?

texts:
Is it possible to use arcpy to open up an editor (like Notepad) to show a simple text file? I'd like to do this in a tool that generates text which is supposed to be copied and pasted into another text file.

If there is a built-in method of showing text without all the other tool output, I'd prefer that.

Having the text already selected would be a bonus.


texts: What are FLOOR and CEILIING in PostGIS designated for?

texts:
Just a short question to lighten up my understanding of some PostGIS-code:

Regarding the creation of user-defined grids I often find the expressions CEILING and FLOOR. Like in this non-representative example:

  SELECT floor(st_xmax(the_polygon)/1500) FROM ...

What are they exactly designated for?

Thanks in advance for lightening me up!

knutella


texts: Why is the map content not displayed when I reopen Print Composer?

texts:
I just updated from QGIS 1.7 to 2.0 on Mac. When I now open a project and create a map in the print composer the content of the map is lost when I open the project again later. The settings of the stored print are still ok, but the map content is not displayed. Instead it just says "Karte wird hier gedruckt" which means "map to be printed here" in English.

Does anybody have an idea what the problem could be? It always worked fine with 1.7.


texts: AxMapControl: How to set map extent in ArcObjects?

texts:
I have a stand-alone application written in C#. I want to be able to save and restore the current map view (extent) that the user sees. I assumed that it would be as simple as saving the active view's extent and then later resetting it, but that doesn't seem to work.

Here is some code:

_myModel.Envelope = _axMapControl.ActiveView.Extent; // save it

...

_axMapControl.Extent = _myModel.Envelope; // doesn't work
_axMapControl.ActiveView.Extent = _myModel.Envelope; // also doesn't work

_axMapControl.Refresh(); // doesn't do anything
_axMapControl.ActiveView.Refresh(); // doesn't do anything
_axMapControl.ActiveView.ContentsChanged();  // doesn't do anything

Note that after the assignments, both _axMapControl.Extent and _axMapControl.ActiveView.Extent values are unchanged.


texts: Looking for high resolution land use (not land cover) data sets for world?

texts:
I am looking for high resolution land use data sets for the world.

This question is NOT to be confused with an earlier question - Equivalent of CORINE land cover data sets for rest of world which refers to land cover.

My question relates to high resolution land use data sets.

http://en.wikipedia.org/wiki/Land_use


texts: I have a shapefile in Wagner VII projection (well, the proj4 string really only is: "+proj=wag7 +datum=WGS84 +units=m"). I'd like to reproject it to WGS84 (EPSG:4326) using QGIS 2.6. However, when I try to do so, the program outputs an error message, telling me that the forward projection was not successful, with the following error: "non-convergent inverse meridional dist". What am I doing wrong?


texts: Reason for a Data Loss While Making Schema Changes on SDE - Oracle DB

texts:
I am not an expert when it comes to Administering Multi-User Geodatabase. However, I refer to Esri's documentation for best practices and efficient workflows before putting on my GIS Admin's hat to perform any Geodatabase Administrative task.

Recently, one my tasks involved making schema changes (adding fields and deleting fields) to a feature class part of feature dataset which is registered as versioned. As per the Esri documentation here, it is required to unregister the dataset as Versioning before making any schema changes to the feature class. Keeping this mind, I followed the steps highlighted below:

Take a Geodatabase and XML Backup of the Production DB
Connected to the Geodatabase Administrative Account
Removed all the locks from the dataset
Connected to the Geodatabase Data-owner Account
Unregister the Dataset as Versioning - ArcGIS Prompted me to compress the edits to the base - I checked option before hitting Continue (I remember checking the mark but have no way to confirm I did)

Added the fields and deleted the fields required to complete the task
Registered the feature dataset as versioned

ArcGIS for Desktop Version - 10.4.1
Oracle Enterprise Database Version - 11g

Later, the data owner reported that was data missing from the user's version. Upon investigating further, the data was missing from DEFAULT Version and User's Version (one year worth). The data was retrieved from the backup but my  still need to know what caused such huge data loss.

My research tell me that the data loss occurred due to not checking the compress option. However, the data editor's version regularly post their edits to DEFAULT which should logically delete only the data (edits worth 3 days) that was not posted to Default before making the above changes.

I am asking for some guidance from experts in the field who can propose any reason for such huge data loss or guide me to the right direction to find a way to look into the System Tables to answer this question.


texts: Write coordinates to new shp file columns with OGR

texts: I'm trying to get the location coordinates written to a shp file after reprojecting it.  I have shp files that are in WGS 84 (EPSG:4326) and I reproject them into UTM 13N (EPSG:27700).  I then want to add two columns and populate them with the UTM coordinates, the new projection.  I know that this is easy to do with QGIS in the attribute table but I have several files and would like to find a way to batch the process and handle several files at once so I am hoping there is an OGR method.


texts:
I'm a novice GISer who needs to do some basic computations on very large (say ~20GB) files. Nothing complicated, merge, intersect, add XY columns kind of thing. It's a necessary step for a project, but it won't be something I do frequently. From what I understand, ArcMap refuses to work with files over 2GB and I think QGIS would crash with these files.

As I can see it, my options are:

Knuckle down and learn GIS python scripting just for this one computation
Break my computation down into 2GB chunks and run it 10 times, hoping none of the intermediates exceed 2GB.
Find software or file formats that can handle large files.

Number 3 would be the easiest if such software exists but I'm not sure it does. If not, do I have any options other than 1 and 2? Both are unpalatable (being new at scripting for 1, tedium and possibility that it doesn't work for 2). Are there other options?


texts: How to parse CRS from shapefile using geotools

texts:
I need to retrieve CRS from a shapefile. I have tried to get it using getCoordinateReferenceSystem() which is provided by org.geotools. My code is as follows.

URL someURL = new URL("https://github.com/basilmohammed/CRS-Extraction-From-ShapeFile/raw/master/ShapeFiles/3857/3857-Point/3857Point.shp");
ShapefileDataStore SHPdataStore = new ShapefileDataStore(someURL);
ContentFeatureSource featureSource = SHPdataStore.getFeatureSource();
ContentFeatureCollection featureCollection = featureSource.getFeatures();
SimpleFeatureIterator iterator = featureCollection.features();
SimpleFeatureType schema = featureSource.getSchema();
CoordinateReferenceSystem sourceCRS = schema.getCoordinateReferenceSystem();
// the coordinates system I want to reproject the data to
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
//Define a MathTransform object
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);

But always I am getting sourceCRS with null as follows.

Any advice?


texts:
I appear to be having problems with re-projecting a bunch of shapefiles that are set to project at WGS_1984_Complex_UTM_Zone_28N coordinate system into Sierra_Leone_1968_UTM_Zone_28N.

I am using the "Project" tool from the Geoprocessing tools to get them into the right projection however after this process is complete the features in my shapefiles are offset. I'm unsure how to correct this, as ArcMap appears to calculate this by default. Does anyone know of a custom transformation that might help with my offsetting issue?


texts: I have raw GPS locations for a track, and need to get maximum speed. Simple methods (calculate speed between 2 points, take max) will give meaningless large numbers due to GPS inaccuracy, the point is jumping around. Can you suggest a good ready-made algorithm to solve this?


texts: How to Hide and Show ArcGIS API ClusterLayer on Demand

texts:
I followed This API Tutorial to create Cluster Points on the map. Now I would like to have a function to hide them. I already tried like?

clusterLayer.hide();

but it didn't work! I also tried to bind the clusterLayer to a GraphicsLayer like below

  [var clayer = new esri.layers.GraphicsLayer();
  clayer.add(clusterLayer);
  clayer.hide();][1]

but still not working
After hours of working I find a sample on google at this Demo doing the job of switching between Points and Cluster but it is supper complicated to me!
Can you please let me know how I can hide and show the cluster layer on demand?


texts: Comparing ranges and distributions of species, looking for correlation and differences between ranges

texts:
As a biology student just beginning with GIS, I wondere if someone could give me an idea of what tools to go about using to compare the range of 3 species that have trophic effects on one another.

I have 3 vector layers, of count data plotted as xy coordinates on a map of the UK, taken from latitude and longitude recorded data.
The first is Water Vole distribution, the second European Rabbit, the third American Mink.

The aim of my assignment is to see whether their ranges coincide (Mink and rabbit supposedly do), or whether they occupy different areas (supposedly, vole and rabbit).

I am looking for a tool that will help with this, but as a GIS beginner, I have explored using point density tool to give me nice density maps in raster form. But am unsure how to get a quantitative output for analysis from this (key for my report).

I am looking to see if there is a correlation/link between rabbits and lynx, and if there is an inverse relationship between voles and lynx. I need a statistics tool that can give me outputs for this.


texts: Generate channel network (rainwater run off) - which algorithm is the best?

texts:
I want to compute DEM for a whole country (cca 79 000 sqr. km) and I need to be hydrologically correct as much is possible. After, for this correct DEM I need to compute drainage lines of surface run-off. I can use 4G model, which is a cloud of discrete points 5x5 m. But this is secondary. Which algorithm should I use to solve this?

On sample data I tested TauDEM, which appears to be a really good tool, but when it's supposed to make raster with drainage lines, the output raster has only one value (whole black). I tried it on ArcGIS 10.2 and 9.3 (on Win 7). Both were the same, first few steps OK, then black raster. I tried to use TauDEM in QGIS, but there was an error:

  "Oooops! The following output layers could not be open"..

I checked paths to toolbox and everything was fine. Next tool that I tried was ArcHydro toolbox, but this tool seems to be very inaccurate.

Does anyone know a good tool to simulate drainage lines? Or why I can't do it with TauDEM?

EDIT:
As a final result I need shapefile of run-off lines. I already have some, but there are many error like crossing of lines, lines aren't snapped to next segment and so on. Now I'm trying to compute them again without these errors.

Example of run-off lines:


texts: Define the "system" in Geographic Information System

texts:
When speaking or writing about GIS all too often the concept in peoples' minds for Geographic Information Systems can be rephrased as "the software for map data". While software is integral to GIS, thinking of spatial software as the GIS  is a narrow view, unnecessarily limiting our conversation and understanding.

I have some ideas about what the "system" in geographic information systems really means, which I'll throw into the punch bowl below, but let us hear yours first.

Secondly, is it too late to bring into public consciousness an expanded truer definition which will arise in answer to the GIS moniker, or do we need to find a new term to encompass it? if the latter, what do you suggest?


texts: New Delhi, India, town level boundary files for GIS

texts:
I'm trying to find boundary files for New Delhi to map religious diversity (data here:
http://www.censusindia.gov.in/DigitalLibrary/MFTableSeries.aspx), and I have downloaded the GADM files and tried some other sources mentioned in earlier threads. My data is at the town level and I am looking for a shapefile at the town level for New Delhi. So far, the files I've found do not match the boundaries of the towns/areas listed in the excel sheet.

Can someone help me find the right file?


texts: Why does Zonal Statistics as Table return an empty table?

texts:
I have a set of polygons with a field of type "SHORT" indicating whether each unit is habitat (1) or not habitat (0).   I also have a raster file in GRID format (32-bit, floating point):  I would like to use the Zonal Statistics as Table tool to calculate two values, the mean value of the raster for habitat areas and the mean value of the raster for non-habitat areas. When I attempt to do so, the tool returns an empty table with no output whatsoever.

I can successfully perform Zonal Statistics as Table using the same polygon layer and 'habitat' field with other rasters in my map document. I have not been able to successfully produce output with this raster and any field in my polygon layer. This leads me to believe there is something about this raster file that is causing this issue, but I cannot figure out what it is. I have read that Zonal Statistics should be able to use a 32-bit, floating point raster as input, so I am wondering if the issue is related to the GRID format or another attribute of this raster I am overlooking.


texts: Create paths layer without duplicating geometrical data

texts:
I have a grid of lines in a geopackage multiline layer (layer_grid). Using the geometrical information contained in this layer I want to create a second layer (layer_paths) that store information (name, description, length, ...) of some paths that can be obtained connecting the lines, but without duplicating geometrical data.

Avoiding the duplication is important since I need to use the obtained database in a mobile application.

How can I perform this task?

I'm using QGIS 3.4.14

EDIT

I'll try to explain with a very simple example.

Lets say that the layer_grid (that contains all the geometric data of the road that i'm considering) is represented by the orange lines. In this layer are stored also the properties of the road, such as for example the type of road (asphalt, dirt, ..). Apparently, every line between two intersections is a different vector in the layer.

In the second layer I want to store the properties of an itinerary between, for example, point A and point B. But I don't want to re-store the geometric data in this level since I need to download these layers (exported in geoJson) from a mobile app.

I hope that I explained the concept.

Clearly, I'll have a greater number of lines and of itineraries and for this reason I need to save as much data as possible.


texts: Identify X Y coordinates of polygons and add them to the attribute table in QGIS?

texts: I have a polygon shapefile which I have imported into QGIS, however it does not display X Y coordinates.  Is there a way that I can identify the coordinate for the centroid of the polygon and display this in the attribute table?


texts:
I am fairly new to all of this and am looking for some guidance. Any information is appreciated!

I am trying to use leaftlet to display shapefiles from an ArcGIS Server using the REST URL. I have seen several examples of this (http://geojason.info/leaflet-vector-layers/demos/arcgis-server/) (http://blog.davebouwman.com/2011/08/04/leaflet-lean-mean-javascript-maps/) and have tried to apply the same formatting to my own code but have not been successful in actually displaying the data on the map. Because I'm not extremely familiar with this stuff I'm not sure what my next step should be to solve the problem.


texts: Problem merging two rasters in QGIS 3.14

texts:
When I try to merge this two or more rasters with QGIS 3.14 this error occurs! Any idea how to solve it?
Versione di QGIS: 3.14.16-Pi
Revisione codice QGIS: df27394552
Versione di Qt: 5.11.2
Versione di GDAL: 3.0.4
Versione di GEOS: 3.8.1-CAPI-1.13.3
Versione di PROJ: Rel. 6.3.2, May 1st, 2020
Algoritmo in esecuzione...
In avvio l'algoritmo 'Fondi (merge)' …
Input parameters:
{ 'DATA_TYPE' : 5, 'EXTRA' : '', 'INPUT' : ['/vsizip/E:/LAVORI/VIA VILLACIDRO FOTOVOLTAICO/FOTOVOLTAICO/DTM/546.zip/50000_546.asc','E:/LAVORI/VIA VILLACIDRO FOTOVOLTAICO/FOTOVOLTAICO/DTM/50000_547.asc'], 'NODATA_INPUT' : None, 'NODATA_OUTPUT' : None, 'OPTIONS' : '', 'OUTPUT' : 'C:/Users/Luca Porru/Desktop/pippo.tif', 'PCT' : False, 'SEPARATE' : False }

GDAL command:
python3 -m gdal_merge -ot Float32 -of GTiff -o "C:/Users/Luca Porru/Desktop/pippo.tif" --optfile "C:/Users/Luca Porru/AppData/Local/Temp/processing_CqXWIW/10f33530bc4c49b0b7bcb98d4c5ac0c0/mergeInputFiles.txt"
GDAL command output:
Traceback (most recent call last):

File "C:\OSGEO4~1\apps\Python37\lib\runpy.py", line 193, in _run_module_as_main
 "__main__", mod_spec)
File "C:\OSGEO4~1\apps\Python37\lib\runpy.py", line 85, in _run_code
 exec(code, run_globals)
File "C:\OSGEO4~1\apps\Python37\Scripts\gdal_merge.py", line 38, in <module>
 from osgeo import gdal
File "C:\OSGEO4~1\apps\Python37\lib\site-packages\osgeo\__init__.py", line 41, in <module>
 _gdal = swig_import_helper()
File "C:\OSGEO4~1\apps\Python37\lib\site-packages\osgeo\__init__.py", line 24, in swig_import_helper
 _mod = imp.load_module('_gdal', fp, pathname, description)
File "C:\OSGEO4~1\apps\Python37\lib\imp.py", line 243, in load_module
 return load_dynamic(name, filename, file)
File "C:\OSGEO4~1\apps\Python37\lib\imp.py", line 343, in load_dynamic
 return _load(spec)
ImportError: DLL load failed: Impossibile trovare la procedura specificata.

Execution completed in 0.23 seconds
Results:
{'OUTPUT': 'C:...../Desktop/pippo.tif'}


texts: Programmatically defining WMS/WFS sources?

texts:
I would like to locate a QGIS project on a shared ressource.

This project makes use of WMS / WFS data.

Is there a way at project opening to programmatically define theses sources (if not existing on the workstation) so that the project gets usable by anybody?

Any hint?


texts: Postgis csv and error loading into mysql

texts:
I want to load spatial data into MySQL using a CSV I have generated from a Postgres table and when loading it into MySQL database using phpmyadmin it complains "Cannot get geometry object from data you send to the GEOMETRY field." How best can I export PostGIS geometry into MySQL using a csv without generating errors?

I have created a MySQL table as

CREATE TABLE cadastre (
  gid int(11) auto_increment NOT NULL,
  geom MultiPolygon NOT NULL,
  lpkey varchar(50),
  SPATIAL INDEX(geom),
  PRIMARY KEY (`gid`)
 ) ENGINE=MyISAM  DEFAULT CHARSET=latin1;


texts: Extracting raster NDVI images for summer months over several years (GEE and LANDSAT)

texts:
I want to extract NDVI raster image that is including only NDVI band over a certain region (noted as point) for summer months over several years. I am working with Landsat 8 and the time frame is June, July, and August between 2013-2019. I have written so far the following code:
batch = require('users/fitoprincipe/geetools:batch')
var cloud_masks = require('users/fitoprincipe/geetools:cloud_masks');
var landsatSR = cloud_masks.landsatSR;

// * Landsat 8 SR Tier 1 2013-2020, B5 and B4
var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterBounds(point)
                  .map(landsatSR());

var ndviFun = function(im){
  var ndvi = im.normalizedDifference(['B5','B4']).rename('NDVI');
  return im.addBands(ndvi);
};

And then map ndviband function over each image within dataset using the following code.
var ndvi = dataset.map(ndviFun).select('NDVI');

And now I am trying to create a time-series loop where each year summer ndvi is reduced by mean and the end result is image collection where one raster image for each summer between 2013 - 2019.
// get a number of year to loop
var nYear = ee.Number(endDate.difference(startDate, 'year')).round();
print(nYear)

// get a list of time strings
var yearList = ee.List.sequence(0, nYear).map(function (n) {
  return startDate.advance(n, 'year').format('YYYMM')
})
print(yearList)

var result = ndvi.map(function(image){
  // map over each month
  var timeSeries = ee.List.sequence(0,nYear).map(function (n){
    //calculate the offset from startDate
    var ini = startDate.advance(n,'year');
    // advance 3 month to take summer
    var end = ini.advance(3,'month');
    //filter and reduce
    var mean_ndvi = ndvi.filterDate(ini, end).reduce(ee.Reducer.mean());
    // get the value and check that it has data
    var val = ee.Number(mean_ndvi.get(['NDVI']));
    val = ee.Number(ee.Algorithms.If(val,val,-999));
    // return mean
    return val;
  });
});
print(result)

I run into an error that is saying "User-defined methods must return a value." Any suggestions?
Here is the code link for a snapshot: https://code.earthengine.google.com/474f6d857870f3698d2cf3db96720083


texts: GEE - reduceRegions for an image collection

texts:
I am relatively new to Google Earth Engine, currently working on finding mean NDVI for command areas of different tube wells in my study area.

What I have done so far is:

Filtered L8 image collection for the months of Jan and Feb for 2014-2019
Imported the feature collection which contains the command areas
Masked out clouds, and mosaiced according to acquired date (which reduced the size of image collection from 377 to 87)

What I am trying to do is, to apply reduceRegions and map it to the image collection, so that I can export a table of 613 (number of features in the Fc, sorted by ID) * 87 (images in the image collection, sorted by date)

When I try to do reduce regions for a single image, the output is a feature collection with 613 features and columns corresponding the the attributes of the features

var reduced = mapNDVI.first().reduceRegions({
    collection:cmdArea ,
    reducer:ee.Reducer.mean(),
    scale: 30
  });

      which gives

featureCollection (613 elements, 46 columns)
type: FeatureCollection
columns: Object (46 properties)
features: List (613 elements)

When I try to map this over the image collection of 87 L8 images,

var reduced = mapNDVI.map(function(image){
  return image.reduceRegions({
    collection:cmdArea ,
    reducer:ee.Reducer.mean(),
    scale: 30
  });
});

      this gives

FeatureCollection (87 elements, 0 columns)
type: FeatureCollection
columns: Object (0 properties)
features: List (87 elements)

Can anyone help me with how to go abou mapping reduce regions over an image collection to obtain a table of (number of features * number of images)?

Please find the code below:
https://code.earthengine.google.com/7dd30af72b3685f673d4d1afa97334c4


texts: Developing visualization on zip code level data

texts:
I'm trying to use the zip code data from the census (any year that works). I have 5-digit zip code data that I need to match up to a map/shape file. The problem is that the zip code tabulation files from the Censsus ZCTA that can be found here http://www.census.gov/cgi-bin/geo/shapefiles2010/layers.cgi have holes in them. I am assuming these are lakes, parks, and wilderness areas. However, every zip code map I've seen do not have these. Is there a resource or method of combining these to create a solid zip code area?

Also, are there any other alternatives to zip code tabulation areas? Google maps will not work for what I am doing unfortunately. :(

EDIT: I am aggregating zip codes based on the data I have. I was hoping that some of these regions would be absorbed by other zip codes so when they are aggregated (based on the data I have) it would show different regions without these holes.


texts: Checking CRS of project layers against reference one

texts:
I have a lot of layers in my project. I need to compare coordinate reference system (CRS) of these layers with another layer used as a reference layer.

I wrote a script to do it, but it didn't work like I thought it would.

mapcanvas  = iface.mapCanvas()
layers = mapcanvas.layers()

for l in layers:
    #print l.name()
    print l.crs().authid() #printando o epsg
    if l.crs().authid() == "ESPG:31983": #here I would check SRC
        print 'ok'
    else:
        pass


texts: How to convert to web mercator projection using an open source tool?

texts: I've got a KML file which I need converted to a web mercator projection. This file will go into Tilemill. What are some good open source tools to convert files to the web mercator projection?


texts: How to generate xml for use with geoxml3

texts:
I have a kml and shp I would like to convert to an xml to use with geoxml3. How would I go about doing such a thing? The xml I generate should look something like this:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2"     xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Document>
<name>Districts</name>


texts:
I have several field in an attribute table that have  value, I am wanting to using the 'Calculate Field' GP tool to select the  values and replace them with zero.

Anyone have any ideas on the Python and/or SQL statement I would use to execute this?

This is the Calculate Field GP

I tried this but it is not working


texts: Create Shapely linestring from string of coordinates in text format

texts:
I have a DataFrame where one column consists of lines stored like so:

'[[-7.61, 42.15],[-6.62, 43.16],[-8.9, 41.0],[-9.2, 40.5]]'

Is there an easy way to read these lines directly into shapely LineStrings so I can treat this as a GeoDataFrame?


texts:
I want to transform about 750 adresses into geocoded points. I´ve tried the plugins MMQGIS (Geocode) and Geopunt4Qgis but none of the things I´ve tried worked. Maybe its because the adresses are german or the table is not formatted correctly before I export it to a CSV file.

the table looks like this:

this is what I type in MMQGIS (the City is at the time a state, there is also no option for housenumbers)

this is the error I get

this is in Geopunt4Qgis (there is an option for housenumbers but I don´t have the postal code)

I don´t necessarily need this to be done in QGIS - I just need a point shapefile so any solutions apart from Qgis are also welcome.

Any ideas?


texts: Bowtie or hourglass polygon validity issue when self-crossing point is not defined

texts:
Using Python2.7 and shapely, let:

import shapely
coords = [(0, 0), (0, 2), (2, 0), (2, 2), (0, 0)]
bowtie = shapely.geometry.Polygon(coords)
bowtie.is_valid

which gives this :

Self-intersection at or near point 1 1
Out[2]: False

So let try the buffer(0) method as describer in the doc (http://toblerity.org/shapely/manual.html):

clean = bowtie.buffer(0)
type(clean)

which gives that:

Out[4]: shapely.geometry.polygon.Polygon

The difference between the example in the doc is that in the present case, the crossing point is not defined in the original coordinates, so it's not resulting in a multipolygon and it's not really consistent with these initial coordinates, as this command:

clean.wkt

shows when it returns:

Out[16]: 'POLYGON ((0 0, 0 2, 1 1, 0 0))'

And the resulting image for a more visual understanding:

With:
Red: original self-intersecting polygon
Teal: cleaned bowtie polygon with buffer(0) method (only one half of the original polygon. Take care of the different scaling effect on iPython...)

Question:
How to easily get the original shape back with a valid geometry so that we can perform further operations like intersections, etc?


texts: Calculated distance doesn't match Google Earth

texts:
I have a google earth route that measures 190.1km in Google Earth.

I then pulled down the .kml file, extracted the points and measured them using Vincenty's formula (ellipsoid earth) and got 190.2.

BUT, if I measure the SAME points along the SAME route to a different distance along that route, I'm way off. Sometimes as much as 2km.

For example, If I want to know the lat and lon of the point closest to 100km, I measure the distance between all the points along the route until the measurement is as close to 100km as it's going to get. When I compare that point to the route in Google Earth, it's around 98km.

Any ideas of why this is the case?


texts: Count number of unique raster pixels in polygon using ArcMap?

texts: I have a polygon shapefile and I want to generate a table of a count for each pixel value in each polygon.  My main goal is to find out what pixel is most abundant in each polygon.  I am using ArcMap 10.


texts:
Please excuse my lingo, I am still learning some aspects of geo spatial and arcpy.

I have a python script that takes a shapefile. The shapefile consists of a polyline. I have traversed the polyline, done some interpolation and retrieved some points. Long story short, I'm trying to generate polygons along the polyline based on intervals of the M values. These polygons have four points (my geometries).

How can I take these geometries, or points, that are in my python script and output them into a shapefile. I know I'll have to pass my points to a function that produces a shapefile, but that is my dilemma. How exaclty do I create the function that creates my shapefile for the geometries?


texts:
I have a shapefile polyline which I would like to add nodes to at set intervals along its length.

Is this possible in ArcGIS 10.2 and how would I go about it?


texts: Intersection between two layers in Atlas Print Composer

texts:
For making an atlas I'd like to show different city districts, each with enclosed water bodies - one district per map (see picture).

I try to show only the part of the matching water geometry which intersects with the the city district or rather which lies only within the district so that the water does not overlap the district's border.

My idea: First step is displaying only water features which have the same district name as the district itself (before the AND-operator), the second step contains the described matter.

But I guess I imposed some wrong conditions within the rule-based symbology of my water layer. Is my expression redundant anyway?

"district"  =  attribute(@atlas_feature, 'district') AND intersection(@atlas_geometry,geometry(get_feature('districts_layer','district', attribute(@atlas_feature, 'district'))))

Does anybody know how to solve this?


texts: Data-driven line thickness in ArcGIS API for JavaScript 4?

texts:
I am trying to render a line featurelayer and would like to set the thickness (width) of lines based on one of the attributes, which is a continuous value between 0-1; for example for value of zero the thickness would be 1, for the value of 1 the thickness would be 10, and other values between 0-1 would be proportionally assigned a thickness (width) between 1 to 10. I am looking at this example from ArcGIS JavaScript API 4.2 in which the size of circles are determined based on an attribute.

Could we do a similar thing for lines' thickness?

In visualvariables what value should I give to the type?


texts:
I am a GIS novice using SQL 2008, with C# and the NetTopologySuite WKTReader type to work with some geometries.

I have some data which I have been able to save into the database but when I read back out the I receive a 24144 error with suggestions to use MakeValid() to convert it to a valid instance.

I don't understand whats wrong with the WKT. I presume its violating an OGC rule but not sure what.

GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (POLYGON ((-2.76 49.38905, -2.7600333333333333
49.535166666666662, -2.5833500000000003 49.5352, -2.4087833333333331
49.535166666666662, -2.4087666666666667 49.38905, -2.56075 49.389099999999992,
-2.5552833333333336 49.437733333333334, -2.5552166666666665 49.455183333333331,
-2.5968 49.455183333333331, -2.5967833333333332 49.437733333333334,
-2.5552833333333336 49.437733333333334, -2.56075 49.389099999999992, -2.76 49.38905))))

I am unable to change this data so is there an alternative to reading this out of SQL without having to change it?

Thanks for looking


texts: Is there an equivalent to arcpy.ImportToolbox() in 9.3.1?

texts: Is there an equivalent to arcpy.ImportToolbox() using ArcGIS 9.3.1? I'm assuming it would use arcgisscripting instead of arcpy but I can't find any documentation.


texts: Iterate on Fire data (FIRMS) in Google Earth Engine

texts:
i loaded the FIRMS data in Google earth engine for the year 2018 and now want to know for a created AOI, when and if there was a fire there. My idea is to create a loop that shows me in a five-day rhythm whether there is a record in the AOI or not.
Does anyone have any idea how to set up something like this?

This is the first part of my code:

//dataset of FIRMS
var FIRMS = ee.ImageCollection('FIRMS')
              .filterBounds(Greece)
             .filterDate('2018-01-01', '2018-12-31');
var T21 = FIRMS.mosaic()
             //.clip(StudyArea)
               .select('T21');


texts: Labeling cells in fishnet using equivalent ("southeast") vertices in grid using ArcGIS Desktop?

texts:
I have a large grid of 0.001-degree fishnet cells. I want to label each cell using the coordinates of its own southeast vertex.

I have the coordinates for all the vertices in the entire grid, but I don't see a way to select the SE vertex for each cell except brute force. I do understand that all the vertices in the entire grid are SE corners of "some" grid cell. My challenge is to associate a single vertex with each cell for label purposes.

I did try identifying the "start", "end", and "midpoint" vertices, but I see that different cells in the grid have each of these vertices at different corners, e.g., the "start" vertex was not in the same corner for all cells, so this didn't work.


texts: crssync.exe error while installing QGIS 3.0 under Windows 7

texts:
I want to use QGIS 3.0 under Windows 7.
Using the OSGeo4W installer I get the following error message:

 The program can't be started because api-ms-crt-runtime-I1-1-0.dll
is missing.

Can I use this installer for Windows 7?


texts: Descending order of VendorOption:sortBy?

texts:
Is there a descending order of VendorOption:SortBy?

What is z-ordering?

  It is a common need to perform road casing against a complex road network, which can have its own z-ordering needs (e.g., over and under passes).

For example: <VendorOption name="sortBy">rank</VendorOption>

the will be drawn as 1,2,3,4,5

I want is it has to be drawn as 5,4,3,2,1


texts: Highlighting selected polygon features without filling colors in QGIS while editing?

texts: When editing in QGIS, how can I highlight the selected polygon features without filling them with colors i.e. just outline the polygon boundary?


texts: display latitude and longitude on googlemaps on click

texts:
I am trying to retrieve latitude and longitude values. When the user clicks on the google map, latitude and longitude values should be displayed in the text boxes.

My code is:

<%@ Page Language="C#" AutoEventWireup="true" CodeFile="sample.aspx.cs" Inherits="sample" %>

<!DOCTYPE html>

    <title>Show Google Map with Latitude and Longitude in asp.net website</title>

    <form id="form1" runat="server">

lat:
lon:

    </form>

But it not displaying anything. Please tell me where the mistake is?


texts: Setting Opacity Dynamically in fill-opacity

texts:
I am trying to use a map I made in Mapbox Studio and data I have locally in js code to change the opacity on different countries on a map.

In my dataLayer, I define the layer as such:

    const dataLayer = {
    'id': 'countries',
    'source': {
      'type': 'vector',
      'url': 'mapbox://url'
    },
    'source-layer': 'layer name',
    'type': 'fill',
    'paint': {
      'fill-color': '#FF0000',
      'fill-opacity':['+', 0, ['number', ['in', 'ADM0_A3_IS'].concat(result).map(x => getOpacity(x)/2) .2]]
    },
    'filter': ['in', 'ADM0_A3_IS'].concat(result)
  };

In the fill opacity param, I call an in command in the mapbox API to get country codes ('USA') and then compare them to a list I have stored in code.

My goal is to use the list of country codes to get the opacity values, which I calculate in JS code. The getOpacity() method converts a country code to a float opacity value [0,1].

Maybe I am not understanding the docs right. How can I map the country code values to opacity, if the opacity data is not stored in Mapbox itself?


texts:
OpenStreetMap lists the malls in Albuquerque, but not the stores
inside these malls.

Using the in-mall map/directory, wouldn't it be easy to add stores
inside the mall?

Is anyone doing this?

Of course, a more fun way would be to visit every store in the mall,
but using the mall directory seems more efficient.


texts: Trying to overlay Ordnance Survey Tiles onto a different OpenLayers TileLayer

texts:
I am trying to add two layers to a Map using OpenLayers 3, one is the information from Bing Maps, and the other is data from the Ordnance Survey data.

I have both maps displaying correctly in separate files, however I am unable to combine the two correctly on one map. When adding the two layers to the same map, it seems the Ordnance Survey layer is positioned at the [0,0] lat long point, rather than over the UK. Is there a way to align the two correctly, so they match each other? The accuracy of the overlay doesn't have to be too exact, as my final plan is to switch between the two layers at different resolutions.

This is my first real experience with GIS, so please forgive any obvious mistakes.

The code I currently have is

    var extent = [0, 0, 800000, 1300000];
    proj4.defs("EPSG:27700", '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717' +
        ' +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs');
    var proj = 'EPSG:27700';

    var resolutions = [2500, 1000, 500, 200, 100, 50, 25, 10, 5, 2, 1];

    var map = new ol.Map({
      //controls: [ new OpenLayers.Control.Navigation(), new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.LayerSwitcher()],
      layers: [
        new ol.layer.Tile({
          source: new ol.source.TileWMS({
            url: 'http://openspace.ordnancesurvey.co.uk/osmapapi/ts',
            params: {
              'VERSION': '1.1.1',
              'LAYERS': '200', // initial value; see view resolution below
              'KEY': 'API_KEY',
              'URL': document.URL
            },
            // needs tilegrid otherwise uses inappropriate global grid
            tileGrid: new ol.tilegrid.TileGrid({
              tileSizes: [200, 200, 200, 200, 200, 200, 200, 200, 200, 250, 250],
              resolutions: resolutions,
              origin: ol.proj.transform([400000, 650000], proj, 'EPSG:3857')
            })
          }),
          extent: extent,
          projection: proj,
          maxResolution: 1000
        }),
        new ol.layer.Tile({
            source: new ol.source.BingMaps({
                imagerySet: 'Road',
                key: 'API_KEY'
            }),
            projection: 'EPSG:3857',
            minResolution: 1000
        })
      ],
      target: 'map',
      view: new ol.View({
        projection: 'EPSG:3857',
        center: ol.proj.transform([400000, 650000], proj, 'EPSG:3857'),
        resolution: 200,
        resolutions: resolutions
      })
    });

    //layers param needs to be changed on change of resolution
    map.getView().on('change:resolution', function(evt) {
      map.getLayers().item(0).getSource().updateParams({LAYERS: evt.target.getResolution()});
    });


texts: View lidR LAScatalog Processing Engine Options settings

texts:
Is it possible to view the settings for LAScatalog processing engine options settings in the current working environment? For example, what is the default setting for opt_chunk_buffer()?
For example when printing no information is given on the default settings
ctg
#> class       : LAScatalog (v1.2 format 1 and 3)
#> extent      : 195000, 330000, 5432000, 5725000 (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-67.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> area        : 31399.17 km²
#> points      : 204.39 billion points
#> density     : 6.5 points/m²
#> num. files  : 31444


texts: Adding more than one contact in FGDC compliant metadata?

texts:
I am new to following the FGDC metadata standards. I have received an SDSFIE compliant geodatabase. My job is to build on this Geodatabase.

Each feature class already has FGDC compliant metadata. As I am adding to feature classes I would like to update the contact information to include both my name/contact and also maintain the name/contract of the previous person who worked on this Geodatabase.

When I go into editing mode via ArcCatalog working on stylesheet FGDC ESRI I have found the "identification" and "contact" tabs. I don't see where I can add an additional Person as a contact.

Is this possible? If so, what am I missing?


texts: In QGIS processing modeller, how to use an input number in an expression?

texts:
I'm making a model where I need to use a certain years in a bunch of expressions. For example:

""era_2" > YEAR1"
"era_2" > YEAR2 AND "era_2" < YEAR1
"era_2" >= YEAR3 AND "ERA_2" <= YEAR1

It's confusing at the moment for other people using the model. So is there a way I can use an input number in the expression?


texts: Selecting features intersecting a selection from another layer in QGIS

texts:
I have two large WFS layers in QGIS (2.16.3): one with polygons and the other with points. I've selected some of the features from the polygon layer and now I'd like to select all the features from the point layer that intersect only with the selected features from the polygon layer. How would I go about doing this?

In the past I had solved the problem by saving the polygon layer's selection as a new shapefile layer. Then I would select the features from the point layer that intersect with the new layer. The problem with this approach is that the size of the shapefiles can end up being very large. I also face the first world problem of having to delete the shapefile once I'm done using it for selecting the features from the point layer.

I had also tried to copy the selected polygons and then pasting them as a temporary scratch layer, but maybe QGIS doesn't like the WFS layer I'm using, or something, because depending on which polygons I've selected, it creates only some of the features for the temporary layer, refuses to create the layer or even crashes QGIS altogether.

I had also tried to use the intersects(), geomintersects() (refFunctions 1.0 plugin) and isselected() (Expressions+ plugin) functions when selecting from the point layer with "Select by expression", but couldn't figure it out.

In another similar question posted on this site it was recommended to use Vector->Research tools->Select by Location with the "Use selected features only" checkbox ticked, but as JaakL had commented, there is no such checkbox for me either. (I would have actually just commented on that question, instead of creating a new question, had I had enough reputation to allow commenting.)


texts: Mapbox - Uploaded GeoTIFF file has no bounds

texts:
I'm trying to add GCP data to an old map and then upload it as a tileset to Mapbox. I've managed to add the data using Georeferencer and the output in qgis looks okay, but when I upload it as a tileset to Mapbox it doesn't show me any preview and when I try to add it as a layer I get an error saying "This tileset isn't visible from your map view".

For reference, "no preview" means:

It gives me these details after the upload:

That seems a little odd to me so I attempted to use the command line tools to find out what's in the TIFF.

If I use gdalinfo on the file (with an argument of -proj4), it outputs what looks like reasonable output for Corner Coordinates:

Upper Left  (  -6.3427176,  53.3848600) (  0d 0' 0.21"W,  0d 0' 1.73"N)
Lower Left  (  -6.3427176,  53.3155539) (  0d 0' 0.21"W,  0d 0' 1.72"N)
Upper Right (  -6.1930269,  53.3848600) (  0d 0' 0.20"W,  0d 0' 1.73"N)
Lower Right (  -6.1930269,  53.3155539) (  0d 0' 0.20"W,  0d 0' 1.72"N)
Center      (  -6.2678722,  53.3502069) (  0d 0' 0.20"W,  0d 0' 1.73"N)

But the proj.4 string is '+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs'.

Are those 0s for lat and long messing things up for me, or am I barking up the wrong tree?

I'm a QGIS noob so step by step instructions might be needed.


texts:
I'd like to create a new shapefile or a new geometry variable that allows me to plot borders around regions in R. I'm using the sf and mapping with tmap. Basically, I'm adding a character vector to an sf object and would like to make the character vector the new/preferred mapping border.

Here is an example of my approach, which doesn't do what I'd like. I can't tell that it does anything.

library(tidyverse)
library(sf)
library(tmap)

## use North Carolina example
nc = st_read(system.file("shape/nc.shp", package="sf"))

nc_new.region <- nc %>% ## add new region variable
    mutate(new.region  = sample(c('A', 'B', 'C'), nrow(.),replace = T))

nc_union <- nc_new.region %>%
    group_by(new.region) %>% # group by the new character vector
    mutate(new_geometry = st_union(geometry)) # union on the geometry variable

# map with tmap package
tm_shape(nc_union)+
    tm_borders()

PS I posted this at SO as I wasn't sure where to put it.


texts: osm2pgsql import planet.osm.bz2 slowness

texts:
I am trying to import whole planet data to my i7 3770k + postgres on SSD + 16gb RAM in order to render my custom tiles.

But the process is painfully slow and multi core is not helping (though you have specified)

osm2pgsql -d planet -U osm --slim -C 8000 --number-process 4 planet-120704.osm.bz2

From htop what I can see is osm2pgsql take 100% of 1 core and postgres takes only 10% of cpu.

and current prompt is:

Processing: Node(208320k 147.0k/s) Way(0k 0.00k/s) Relation(0 0.00/s)

Q: Is my computer too slow? Any method to speed up (without complexing too much hardware)?


texts: Fill Features of attribute table by using another feature in QGIS

texts:
I would like to fill up an attribute table in QGIS for the feature 'Color' based on another feature in the attribute table.
At the moment it looks like this for the features 'Direction' and 'Color':
Direction   Color
north       none
south       none
east        none
west        none
I would like to fill the Color feature with RGB color codes to directly plot them by using my code. At the moment I only managed to add a new Field for the feature Color (the other features are already implemented). For example if there is North, the color code '0,165,0,255' shall be inserted into the other feature.
# b is the shapefile applied

   from PyQt4.QtCore import *
        b.dataProvider().addAttributes([ QgsField("Farben", QVariant.String) ])
        b.updateFields()


texts: How do I centre my map after clicking on a marker in carto js? I'm using createLayer but not sure where my listener should be applied (i.e. to map object created using Leaflet or within createLayer call).


texts: Add (XY) Data in ArcMap is not finding folder connection?

texts:
I'm trying to add XY data to ArcMap but it is working.
I tried several times to connect to the folder but it is not being connected and the windows appears empty afterwards

Any suggestions to solve this problem?


texts: Tiling image to be used with Google Maps?

texts:
I have a 2880 x 2880 image that I would like to tile so that it replace the entire usual world map.

So far I got this working

gdal2tiles.py -p raster -z 0-6 5KNewMap.png abc3/

but for the zoom level 0 the tile is not centred and there is a kind of white border around.

What is the best and easiest way to tile my picture so that It replace the entiere world map?

Also I'm getting really confused with the difference between "raster" and "Mercator" ? If I use "Mercator" it asks me for some points but I don't which one to set so that it covers all the world.


texts:
how to download a file from geoserver gml from a url using javascript?

I'm doing a web application with javascript, but I could not download files from geoserver, the url given entry does not bring any result.

enter the url using this instruction:

new Ext.form.TextField({
                disabled: false,
                        fieldLabel: "URL GML",
                        value:'',
            id:"Ugml"

to save a link:

var Ugml = Ext.getCmp('Ugml').getValue()

to fetch the data:

xmlhttp = new XMLHttpRequest();
xmlhttp.open("GET", Ext.getCmp('Ugml').getValue(),true, User="admin",Password="geoserver");
         xmlhttp.send(null);
         alert(xmlhttp.responseXML.xml);

any ideas

Thanks for your answers


texts: saveACopy in ArcPy, failed to save document

texts:
I have a script that saves a copy of the current .mxd to a location given to the script by a user. However I keep getting an error that says it failed to save the document.

I have writing rights in the folder and I have tried different locations as well, which doesn't seem to work either.

I don't seem to be able to find what I am missing, would anyone be able to tell me where I am going wrong? I am using ArcGIS 10.3.1.

Old code:

import arcpy

import arcpy.mapping as map

Project_number = arcpy.GetParameterAsText(0)
Folder_location = arcpy.GetParameterAsText(0)

mxd = arcpy.mapping.MapDocument("CURRENT")

mxd.saveACopy(r"{}\{}\{}.mxd".format(Folder_location, Project_number, Project_number))

del mxd

I have changed my code a bit with the help of the comments underneath. It works fine when I put a string for "Project_number" and "Folder_location". However when I set it as a parameter, it fails to save.

The parameters would be this:

Folder_location: U:\Example

Project_number: test

So, when I put these in the script as strings, they work. But when I set these as parameters to be chosen when the script starts (with the exact same data) it fails.

Any ideas?

New code:

import arcpy

import arcpy.mapping as map

import os

Project_number = arcpy.GetParameterAsText(0)
Folder_location = arcpy.GetParameterAsText(0)
Corrected = os.path.normcase(Folder_location)
MXD = ".mxd"
MXDName = Project_number + MXD
Pathname = os.path.join(Corrected, Project_number, MXDName)
CPath = os.path.normcase(Pathname)

mxd = arcpy.mapping.MapDocument("CURRENT")

mxd.saveACopy(CPath)

del mxd


texts: How do I serve OSM map tiles (mbtiles) from a Raspberry PI 2?

texts:
My main processing machine is a Ubuntu box running 14.04 with PostGIS and TileMill. I generated my tiles and the resulting mbtiles export is about 250MB (megs, not gigs). The PI runs python, has a storage card big enough to keep the mbtiles export, and wifi. At most, I might have 2-3 devices accessing the endpoint requesting tiles. Typically, only one device will be requesting tiles.

Keeping in mind, I prefer simplicity with minimal hardware requirements over a feature-rich server, what are my options for serving up the tiles from the mbtiles with a PI 2?


texts:
I am having a difficult time figuring this out.

I installed the QGIS 3.10 with the all-in-one package, and there was no error during installation.

However, when I open QGIS, I don't see the processing toolbox, and processing does not exist in plugins. By the way, Python console does not exist in the Plugins menu.

How do I solve it?


texts: Printing to Scale only works for PDF format

texts:
I often need to print out 1:1000 scale drawings to use as a basemap for landscape surveying in the field. We draw directly on to the map so it has to be perfect.

I've found that when I export to pdf it prints out perfectly to scale, but if I try any other image format (png, bmp, tiff, etc.) it prints wrongly (i.e. I measure the scale bar and it's not the right length). I've looked around but I can't see any reason for this in the help.

I can't print directly from the project as these files needs to be emailed to other colleagues for them to print for themselves. I've double checked that I'm not using 'Fit to Page' or any similar other printing related things. So, my question is - is PDF the only possible format for this officially, or are there other formats that work properly?


texts: Enlargement in one dimension to bigger polygon

texts:
I have areas (polygon) in the streets (polygons).

I want to enlagre the parallel segment to an other polygon to touch it. And so creating this new polygons.

There is a lot of new alogithms in field calculator. Maybe processing alorithms could be a solution or Python or GRASS. I don't know which way to choose.

At this time, i try this way :

Cut areas by street polygons
Join attributes by location to grab street polygon id in the area
Polygon to lines
Explode lines
Maybe ** IF segment is < 0-45 degrees THEN segment is extended to road polygons ( exploded too ?)**

Polygons to lines


texts: GeoJSON features are added to null island when dynamically loaded in Openlayers3

texts:
I'm trying to add features from a GeoJSON response in openlayers3.
Using the default loading method works great:

layer.setSource(new ol.source.Vector({url:url, format:new ol.format.GeoJSON()});

But I need error handling, So I use the response from a standard jquery loader like so:

var features = new ol.format.GeoJSON().readFeatures(response);
layer.setSource(new ol.source.Vector({features:features}));

All the features appear bunched on null island in a sort of a mini-map:

And a closer look reveals the features:

edit:Specifying the data/feature projection doesn't work either.

EDIT:
I finally fixed it, Here's the code:

var reader = new ol.format.GeoJSON({
    defaultDataProjection:'EPSG:3857',
    Projection:'EPSG:3857'
});
layer.getSource().addFeatures(reader.readFeatures(res,{
    dataProjection:'EPSG:4326',
    featureProjection:'EPSG:3857'
}));

What an absolute nightmare that was.


texts: Polygons Shrinking and Moving as if CRS Was Suddenly Out of Sync with Project

texts: I long ago set the CRS for a linguistic geography project I'm working on involving Oklahoma data (EPSG:2836 - NAD83(HARN) / Oklahoma North), and I have been closely monitoring that the polygon shapefiles I have been creating are set to the same. I've been working on the project for several days now, and all seemed fine today when I started working. But after hours of work this evening, I saved my project again only to discover that all of the polygons that I made before today have now moved hundreds of miles away to the southwest (while still maintaining their relative locations with respect to one another) and shrunk to ridiculously small size. For instance, a nearly 220 km wide polygon originally in northeastern Oklahoma is NOW ONLY 3 METERS WIDE and located in Southeastern New Mexico! All of the shifted polygons still have the same CRS settings as the others, have same relative position with respect to one another, and their attribute tables are intact. I reiterate that every polygon I've made so far in this project is its own separate shapefile, most of which overlap with one or more others. Oh, yes, I freely admit that I'm very new to QGIS, and GIS in general, only having started to teach myself the basics this summer. But I would GREATLY appreciate figuring out what I did wrong and how can I fix it (beyond just scrapping those polygons and starting over again). I don't know of anyone else that uses the program, so there's no one I can turn to.


texts: Extracting by Attribute Using ModelBuilder with User Input?

texts:
Using ModelBuilder, I want to create a tool that will let a user choose the value or values that will be used to extract a shapefile.

I have a parcel layer with a "Parcel_Type" field. I want the user to be able to choose the type of parcel to extract. So if the user wants to extract the parcels with a type "EL" then he would only have to enter "EL" as a variable. If he wants type "EL" and "CDD" he needs to enter them as multi-variables. That would be the only human interaction with the model.

So how can i use those variables in my model?


texts: Ways to access WMS map server from Qt QML "Map" component

texts:
I have an application which uses Qt library (version 5.8). Said application needs to draw a map, provided by a (local network) map server. The server gives only WMS maps (and is maintained by a third party, so "WMS only" clause is a given).

Qt provides a standard QML cartographic component, Map. It has several out-of-the-box plugins (which basically are responsible for actual tile requests, local tile cache maintenance and some other stuff). Surprisingly, the documentation (see paragraph "Plugin References and Parameters") doesn't tell anything about WMS support (it uses OpenStreetMap as "default" option). Even more surprisingly, information on the subject seems to be scarce. Here are the options I found:

There is ESRI QML class for displaying WMS layers, which requires installing the whole ArcGIS SDK, which looks heavy;
One can write a QML plugin of their own. But even though all I need is to change the request URLs, there seems to be no easy way to reuse existing logic of, say, OSM plugin for other plugin responsibilities (and that's a substantial amount of code);
A web browser component can be used in QML; it can open a local HTML page with OpenLayers or the like. That's significant overhead and in my particular case will cause problems with proxy settings (the browser component is already used in other part of the application).

My question is: are there any other known ways to use WMS server with QML components? Cartography isn't the primary focus of the application, so it doesn't require complex features like geocoding or routing.


texts: Setting output filenames from QGIS modeler?

texts:
I have built model which iterate through features in input layer.

How I can set output filenames in the modeler using input rows?

Something like current_field_value.shp


texts: How to change the default basemap in a Web AppBuilder app via code?

texts:
I created an app in the Web AppBuilder v2.0, running at http://my-PC:3344/webappbuilder, and included the standard BaseMap Widget.

I downloaded the resulting app and deployed it using IIS at C:\inetpub\wwwroot\MyApp. I've subsequently made some changes by editing the JS files directly in an IDE (not via the Web AppBuilder UI) and it's all working fine.

Now I'd like to change the default basemap which is shown when the app opens.

\MyApp\config.json lists the BaseMapGallery widget
\MyApp\configs\BasemapGallery\config_BaseMap Gallery.json lists the
available basemaps
neither file seems to configure the default basemap

Which files should I edit to change the default basemap?

Edit: I am using the LocalLayer widget to specify the map layers manually


texts: Get single part polygon by their intersection in PostGIS

texts:
I have a table in PostGIS that has polygons overlapping each other. (see below an example of 4 overlapping polygons )

Someone could give me an idea, what PostGIS function I can use that give me single part polygons based on their intersection, like the following:

PS: Actually what I want in the end is the individual polygons with no overlap, but the count of polygons overlapping it (counting itself). For example, in the second figure, the count of the polygon 2, 8, 12, 13 will be, respectively, 1, 2, 3, 4


texts:
I have a points dataset (X and Y values) with different values (from 0 to 5). I wish to obtain a map like that:

I need the points with a value of 5 to be represented with the darkest color, the point with a value of 4 with a lighter color, etc.

How can I do that with QGIS?


texts: Is it possible to extrude a vector layer with QGIS 2.16 Globe plugin?

texts:
With the new version of Qgis ( 2.16) it is stated that the extension was rewritten globe .

Is it possible to extrude a vector layer with a height parameter and if it is possible to have the procedure?


texts:
I am using the OREGONSTATE_PRISM_Norm81m data set on Google Earth Engine and I am trying to clip the data to only Washington State. I uploaded a shapefile from the Washington Department of Natural Resource GIS Open Data. This is the code I currently have, I only have added the shapefile and changed the map center from the original data.

var dataset = ee.ImageCollection('OREGONSTATE/PRISM/AN81d')
                  .filter(ee.Filter.date('2018-08-01', '2018-08-15'));
var precipitation = dataset.select('ppt');
var precipitationVis = {
  min: 0.0,
  max: 50.0,
  palette: ['red', 'yellow', 'green', 'cyan', 'purple'],
};

var state = table.filter(ee.Filter.eq('JURISDIC_3', 'Washington State'));
var rain = precipitation.clip(state);

Map.setCenter(-120.332, 47.448, 7);
Map.addLayer(precipitation, precipitationVis, 'Precipitation');
Map.addLayer(table);

The error message states that "precipitation.clip is not a function" and I am very confused because I have used a similar code several times before.


texts: I am working on a field survey and a client requirement is to differentiate various norths, like magnetic north, true north and grid north. My question is what is the methodology of calculating north and how can we differentiate various type of north as above mentioned.


texts: ArcGIS - Percentage in Census tract to HUC

texts:
I’m trying to figure out the percentile of demographic value (minority %, poverty %, etc.) per watershed at HUC12 level. I was hoping that I can somehow convert the values from census tracts into HUCs.
For demographic value dataset, I am using CDC’s Social vulnerability index, which has population estimate, percentage, and percentile for each of the demographic value by census tract.
Initially, I thought I could use the “Summarize Within” function on ArcGIS by using Area-Average calculations; however, this doesn’t make sense for what I’m trying to solve since this function calculates the average per census tract within the HUC rather than the total within the HUC.

Red outline is the watershed boundary. Blue with black outline is the census tract.
Does anyone know of an approach to use to calculate the percentage of demographic value in each HUC using census tract level data? I’ve thought of intersecting the census tracts into the HUCs and merging them back to find the weighted average, but I’m wondering if there’s a more efficient way to do this.


texts: How do I create multiple polylines from single polyline that consists of disjunct sections?

texts: My goal was to create a network of polylines in NHD plus based on start and end points associated with sampling locations. Dissolving by stream name was my first step and has worked just fine for named streams, as I just split them accordingly in an edit session afterward. However, after dissolving, all of the unnamed sections are represented by one entry in the attribute table of the new dissolve shapefile. This was not a big deal until I realized that some sites were on unnamed sections of the NHD network.  When I select an unnamed segment in an edit session and try  to split it, an unnamed section 20 miles away is also highlighted. Is there a way to extract single pieces from a polyline that has numerous disjunct sections? Or do I have to retreat back to the pre-dissolved NHD layer and sort out the unnamed sections from there? I suppose I could go back to that original layer and provide names for all unnamed sections of interest, then dissolve by name again but that seems extremely tedious. Thanks,


texts: Why can't I import my shp file into a postgis database with SRID 102631 using shp2pgsql?

texts:
I have been trying to figure out the reason shp2pgsql isn't working when I try to import a shapefile with SRID=102631. I've had no problems with other SRIDs.  For all shapefiles in 102631, I get this error: Failed in pgui_exec(): ERROR:  AddGeometryColumn() - invalid SRID

Here is a link to one file, for example, that won't import.

docks.shp.zip


texts: Conditional Statements in Arc Map 10 (field Calculator)

texts:
I'm attempting to calculate the values of a field based on IF THEN statements in ArcMap 10 (in field calculator) but I haven't succeeded yet on that.

What I want to add is the following expression:

IF [FIELDNAME1] = [FIELDNAME2] then

[FIELDNMCH] = "no change" else

[FIELDNMCH] = "change"

END IF

The values in the fieldnames are texts.

Any help on that?


texts: Data masking and plotting issues xarray/rioxarray | geojson -> geoviews

texts:
I'm loading a NetCDF file (ERA5 surface temperatures), then I mask the data using country shapes, and finally generate a plot using geoviews. For the unmasked data everything works as intended. As soon as I clip using rioxarray, half the globe goes missing. In detail:
# Read all relevant data:
data0 = xr.open_dataset("./data/COPERNICUS/adaptor.mars.internal-1621870453.0792851-5880-10-bd60c7b4-8e84-4d2a-acfb-d2cc442d24e1.nc")
sf = geopandas.read_file("./data/NATURAL_EARTH/ne_50m_admin_0_countries.shp")

# Filter and set CRS:
dataslice = data0.load().where((data0.time.dt.strftime("%Y") == "2020") & (data0.expver == 1)
                                  , drop=True).sum("expver").resample({"time": "Y"}).mean()
dataslice.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
dataslice.rio.write_crs("EPSG:4326", inplace=True)

Plotting this with geoviews looks ok:
dataset = gv.Dataset(dataslice, ['longitude', 'latitude', 'time'], 't2m')
images = dataset_clipped.to(gv.Image)
images.opts(cmap='viridis', colorbar=True, width=900, height=600,
            projection=crs.PlateCarree())

Now I want to clip the data using the shapefile and plot the clipped data:
clipped = dataslice.rio.clip(sf.geometry, drop=False)
dataset_clipped = gv.Dataset(clipped, ['longitude', 'latitude', 'time'], 't2m')
images_clipped = dataset_clipped.to(gv.Image)
images_clipped.opts(cmap='viridis', colorbar=True, width=900, height=600,
                    projection=crs.PlateCarree())

While the data is clipped by the landmass as intended, all of the Western hemisphere is gone as well.

Remarks:

This does not depend on the projection parameter of geoviews.
Both the gridded data and the shapefile have CRS EPSG:4326.
The Shapefile is complete:

Questions:

What is going wrong?
Why do we need to set the CRS anyway? Both the xarray and the geojsons are just tuples (lat,lon,value) and (lat,lon), respectively. Yet dataslice.rio.clip raises a MissingCRS exception if dataslice.rio.write_crs is skipped.


texts: How to convert SRID's on the fly and via command line, on INSERTION's query in PostGIS?

texts:
We are working on a project that used shapefiles in order to build our first PostGIS database, through a tool called “PostGIS Shapefile Import/Export Manager”.

The problem is that now we need to create new tables onto that same PostGIS DB, through executing SQL commands line, as detailing bellow.

0 – VERSIONS:
a – PostGIS: 2.0.7 r13406
b – Geoserver: 2.12.1

1 – CREATING TABLE QUERY:

CREATE TABLE tabela (gid serial, nome varchar, geom geometry(Point, 32724));

2 – INSERTING A POINT QUERY:

INSERT INTO tabela (nome, geom)
VALUES('Teste Renan', ST_GeomFromText('POINT(-40.388527 -20.339450)', 32724));

By inserting that point via command line, Geoserver cannot correctly calculate “Lat/Long Bounding Box”, when I press the “Compute from native bounds” button, at the “New Layer” screen.
It give strange values, denoting some problem.

 P.S.: I have properly set up EPSG 32724 (WGS84 24S).

By doing this process, when I try to preview that specific layer, nothing appear.

Another strange thing that I have noticed is different values for GEOM field, for the “almost” same place:
a - When created by “PostGIS Shapefile Import/Export Manager”:

0101000000DCD8FFFF83BD1541B1FE00804F915D41

or using st_asewkt command:

POINT(356192.999999417 7750974.00006072)

b – When created by command line:

0101000020D47F0000774EB340BB3144C090A0F831E65634C0

or using st_asewkt command:

SRID=32724;POINT(-40.388527 -20.33945)

Any thoughts on how to solve this problem or how to properly create the query for add point to TABELA table?


texts: ArcMap 10.7 - Using Python to do proper rounding instead of truncation of a raster. Can't figure out why an error is thrown

texts:
My code is throwing the following error:

  TypeError: must be str, not float.

It is supposed to use the Int function in Spatial Analyst to do proper rounding instead of truncation. How can I identify my error? I am new to Python and learning as I go. I show below where the error is thrown.

# Name: intrndthm.py
# Description: Converts each cell value of a raster to an integer by proper rounding instead of truncation - will not work for negative values
# Requirements: Spatial Analyst Extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "h:/test"

# Set local variables
inRaster = "demtest1.tif"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute Int
outInt = Int(inRaster + 0.5) - error thorwn here. If I use it in the raster calculator it works. If I Use outInt = Int(inRaster) it works but it truncates.

# Save the output
outInt.save("h:/test/output/demtest1rnd.tif")


texts: I work for a small civil engineering firm in Florida and part of our due diligence on our sites was checking the soil types. Previously, I was using FGDL (a Florida based GIS data clearinghouse) where it was just easy to download the soils data by county and bring it into our maps that way. Well, I went to download some new counties recently and discovered they no longer carry NRCS (National Resources Conservation Service) soils data. I tried going directly to NRCS and download them myself, only to discover that someone at FGDL was linking the spatial data to the tabular data for me and I have zero idea how to do any of that. I've spent most of my day today after lunch trying to figure it out and now I'm desperate for some help. I don't have access to Microsoft Access, which it looks like one might need to connect the data. So I am not sure what to do next. I am a team of 1 at my office, so no one would be able to offer me answers on site. So, How do I get workable Soil Data into ArcMap?


texts:
I have been working on openlayers 4 and specifically my question is regarding limiting the extent the user sees, let's say 2 km by 2km is the viewing area, this configuration was done using the ol.View class's extent attribute.

In that case for the initial zoom level, I'm able to see the entire area of 4kms by panning across, when I zoom into the next level, then the area I can pan across is less than 4kms, how can I maintain the original extent, even on other zoom levels while panning?

Here is the code :

var matrixIds = new Array(14);
for (var i = 0; i < 14; ++i) {
    matrixIds[i] = "EPSG:27700:" + i;
}

var resolutionsArray = [896.0, 448.0, 224.0, 112.0,
    56.0, 28.0, 14.0, 7.0, 3.5, 1.75,
    0.875, 0.4375, 0.21875, 0.109375
];

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.WMTS({
                name: "OS Maps WMTS",
                url: "<layer_url>",
                layer: "Road 27700",
                matrixSet: "EPSG:27700",
                format: 'image/png',
                projection: "EPSG:27700",
                tileGrid: new ol.tilegrid.WMTS({
                    origin: [-238375.0, 1376256.0],
                    resolutions: resolutionsArray,
                    matrixIds: matrixIds
                }),
                style: '',
                crossOrigin: "Anonymous"
            })
        })
    ],
    view: new ol.View({
        center: [529449, 189252],
        zoom: 13,
        maxZoom: 16,
        minZoom: 13,
        extent: [529449 - 1000, 189252 - 1000, 529449 + 1000, 189252 + 1000]
    })
});


texts: Where can I obtain data on flood sites, land sewage systems, and property value?

texts: I am a college student trying to obtain this data for a class project. If anyone can help that would be great. Basically, my hypothesis for this project is that land sewage systems that have not been recently upgraded lead to flooding, thus, lead to a decrease in property value. I would like to obtain data in Chicago, but any other suggestions are welcome.


texts: layer.SetAttributeFilter() does not return correct number of features

texts:
I have been using GDAL/OGR's layer.SetAttributeFilter() in a Python code to select a subset of features in a layer and save to an output shapefile. However, I am not getting the result I expect and cannot figure out why.

Here I outline my process with an example and show why the result is unexpected:

# First, gather a list of features that meet a specific criteria by their FID
keepFIDs = []
# iterate through features in a layer (polygon .shp):
for feature in layer:
    # if feature meets criteria:
        keepFIDs.append(feature.GetFID())

# Then, with keepFIDs, get the SQL query string:
query = "FID IN {}".format(tuple(keepFIDs))
print len(keepFIDs) # = 3991

print layer.GetFeatureCount() # = 7347 total features before applying filter

# Lastly, apply filter to layer and save to output shp
layer.SetAttributeFilter(query)
print layer.GetFeatureCount() # = 551. Why so few?! Should be 3991
print len(query.split(',')) # = 3991, for extra measure to confirm the SQL query had the correct number

# Save layer to shapefile

SetAttributeFilter() should be returning 3991 features, right? I mean, we know that the FIDs in the query exist in the layer I'm filtering because I just got them when iterating through features in the layer. I know from my print statements that the number of FIDs in the "FID IN " is 3991, not 551.

So what could possibly be going on?

EDIT. Here is another example to show bizarre behavior. This is a different shapefile --> layer shown here, with a lot more features:

# With keepFIDs, get the SQL query string:
#query = "FID IN {}".format(tuple(keepFIDs)) - not doing this here. see note below*
print len(keepFIDs) # = 7492
# *note that in this case, query is "FID IN (firstGroupOfFIDs) OR FID IN (secondGroupOfFids)"
# because of parsing issues with more than ~5000 features
# [http://www.samuelbosch.com/2016/08/workaround-for-sql-expression-parsing.html]
# but the output should be selecting 7492 in total:

print len(query.split('OR')[0].split(',')) # --> 4800
print len(query.split('OR')[1].split(',')) # --> 2692

print layer.GetFeatureCount() # = 89550 total features before applying filter

# Now apply filter to layer
layer.SetAttributeFilter(query)
print layer.GetFeatureCount() # = 9382. Why so many?! Should be 7492


texts: How do I step into the C++ functions from an external Python app I'm debugging?

texts:
I just wanted to see how to step into the C++ class QgsVectorLayer()'s methods from the python call to QgsVectorLayer(self.uri.uri(),'foo','postgres').  I'm trying to debug and it trying to step into the function doesn't work.

Essentially what i'm asking is how to debug an external Python app by stepping into the underlying SIP/C++ code?  Are there instructions anywhere?

Or even simply remote debugging the python app by attaching, using visual studio, on a specific QGIS binary...?


texts: ArcGIS Online Not Accepting CSS for Transparency

texts:
I am attempting to make the background of an element (timecontainer) transparent by using ArcGIS Online's custom CSS option.

The desired effect is seen below (note the semi transparent background):

Using inspector in Firefox, I was able to find the element (timecontainer)that needed to be adjusted (below).

The background-color should be changed to rgba(92, 92, 92, 0.41), as seen below.

Using a CSS editor, I drafted the following:

#timeContainer {
  width: 750px;
  height: 200px;
  background-image: url("http://imagebackground");
}

.timeInfo .bg .fg .window-bottom-left {
  background-color: rgba(92, 92, 92, 0.41);
}

I then copied and pasted that CSS into the Time App's custom CSS window.

But the background remained a solid, grey color.

Is there something I am missing?


texts: Reproduce raster with same extent, and same number of columns and rows?

texts:
I have two raster A and B. A and B have the same grid size but different extent.

How can I make A have the same number of columns and rows, as well as the spatial extent of B?


texts: Is there a way to produce a line on a raster, which direction is based on cell values?

texts:
I am looking at a raster layer where each cell has one of 8 values, indicating where flow from the cell is directed. Is there a way to select any given cell, and have a line drawn which will follow the direction of cells as it reaches them.
This can be thought about like a hydrology flow map essentially, using flow direction to draw a path along cells until we reach a pit or output cell. Is there any code out there which may help? I'm sure this must have been done before somewhere.