texts: Has anyone successfully implemented a web service that sits on top of ArcSDE (SQL Server 2008)?

texts:
Here's the situation. We have planned point locations that lie within predefined polygonal project areas. Our IT department currently has these polygonal project areas in SQL Server 2008 as geography objects. When a user enters a point location in a interface, a query gets issued to the project areas table to return the project area that the new planned point location lies within. Works like a champ. However, these project area boundaries change of course in our GIS and then IT doesn't get the changes and we get data busts. IT has talked of writing a web service that would sit on top of our ArcSDE database (again, SQL Server 2008) so they could just make calls to the fresh data we always have in SDE. This of course would require us to store those datasets in SQL Server geom/geog datatypes.

My question is this: Has anyone successfully implemented anything like this? If so, could you share your thoughts. If you have tried and failed, I'd like to hear about that as well. Thanks.


texts: Converting GeoJSON to projected coordinates in Leaflet?

texts:
We are migrating existing angular app from OL to Leaflet. In OL we were able to get the geoJSON in UTM with EPSG:3857 on draw finish However when we complete a draw in leaflet and try to get the geometry and geoJSON it gives the data in LatLong format.

This geometry and geoJSON is getting Saved in DB and considers LatLong as UTM due to string format.

this.map.on(L.Draw.Event.CREATED, function (event)
{
  let layer = event.layer;
  let drawnGeometry = event.layer.toGeoJSON().geometry
  console.log(drawnGeometry);
  let data = drawnGeometry.coordinates;
  self.Coordinates = data.join(); // coordinates are in LatLong format.
  self.cordiData = JSON.stringify(event.layer.toGeoJSON());
});

Please suggest if there is a way out in leaflet to get the Geometry and geoJSON in UTM format.


texts: OpenGeo Suite Community Edition 4.0 Dependency Problems on Ubuntu 12.04

texts:
The OpenGeo Suite 4.0 has some issues with dependencies specifically with respect to GDAL and PostGres libraries.

The issue arises when one tries to install libgdal-dev, which is required, among other thingds, to get ECW reading working.

I tried installing the OpenGeo Suite 4.0 on a fresh install of Ubuntu Server 12.04, and also tried an upgrade from the 3.2 version of the suite. Here's a pastebin of the error from Apt: http://pastebin.com/3BqZVvCn

After an upgrade, ECW's can't be used, event after installing the extensions into the GeoServer's WEB-INF/lib folder. ECW appears as a format, ECW is listed by GDAL as available, creating an ECW store works, but the store is always empty... Valid stores that worked in 3.2 dont' work on 4.0.

A secondary issue is that some file permissions seem to be broken on upgrading from 3.2 to 4.0, these are exhibited by failure to load some aspects of the GeoServer interface, such as logging.

Does anyone else have these issues, and if so, are there any known work-arounds? It seems to me that a resolution would be for Boundless to add in the missing packages to their repository, but that may be naive.


texts: GeoSpatial queries and indexes in memory

texts:
I am currently using Spark/Scala to process some data, and I am at the point that I need to run ~3M of geospatial queries. One common solution would be using MongoDB or PostGIS (or Magellan for Spark which I already know it has not what I am looking for at the moment).

The queries are the classical Proximity queries for GeoJSON/WGS84 to find out, given a point "A", what are the 100 nearest points to "A" in a radius of 1 Mile within a set "B" of similar points.

Is there any library that helps you to do this sort of queries in memory and a library that can create (always in memory) indexes (2DSphere/Geohash?) for the set "B"?


texts:
I made a map with a horizontal legend to show 10 categories in QGIS 2.18.6 (see below), and here comes some questions. Hopefully I can have some suggestions here:

First, I indeed create the legend with 5 columns, but the number of rows is not 2. The right two columns always have three and one categories, respectively. How can I make the legend like a 2X5 matrix?

Second, those color blocks in the legend are arranged vertically. I'd like to have them arranged horizontally; hence, the legend can have all red blocks on the top, and all blue blocks on the bottom. However, I cannot find any setting in the print composer to change the way of arrangement.


texts: Weird XY-resolution for map in ArcObjects

texts:
I have an MXD file with 240 feature-layers, all having the exact same spatial-reference and spatial resolution. I checked this using the following code:

var sortedLayers = allLayers.OfType<IGeoFeatureLayer>().Select(x =>
{
    var featureClass = (IGeoDataset) x.FeatureClass;
    return new
    {
        LayerName = x.Name,
        FeatureClass = featureClass,
        Resolution = ((ISpatialReferenceResolution) featureClass.SpatialReference).XYResolution[false]
    };
}).OrderBy(x => x.Resolution);

This returns 0.00001 for every layer within my map. However when I query the XYResolution for my map I get a really weird number:

Console.WriteLine(((ISpatialReferenceResolution)myMap.SpatialReference).XYResolution[false]);

This code returns 0.0000000024426994295885463. First I thought it is a rounding-issue. But the returned resolution is not even near the expected resolution of the layers. In fact the nearest decimal that can be represented as binary is 0.00000953674, which is 0.0000000000000000101 (=2^-17 + 2^-19). As far as I know the maps resolution is picked from the layers within that map. However there´s no direct way to query a layer's spatial reference and thus its resolution, as the SpatialReference-property on ILayer is write-only.

Has ayone an idea why the map's resolution is different from my layer's resolution?

The spatial reference of all my layers - or more accurate of their feature-classes - is DE_DHDN_3GK2, which is a German
Cartesian coordinate-system on a Bessel elipsoid having metric units. In fact pretty similar to UTM32, however not in tranversal projection. From ArcCatalog:

DE_DHDN_3GK2
Berechtigung: Benutzerdefiniert

Projection: Gauss_Kruger
False_Easting: 2500000,0
False_Northing: 0,0
Central_Meridian: 6,0
Scale_Factor: 1,0
Latitude_Of_Origin: 0,0
Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_Deutsches_Hauptdreiecksnetz
Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_Deutsches_Hauptdreiecksnetz
  Spheroid: Bessel_1841
    Semimajor Axis: 6377397,155
    Semiminor Axis: 6356078,962818189
    Inverse Flattening: 299,1528128


texts: Geoprocessing>Merge and arcpy.Merge_management() producing different results?

texts:
In ArcMap 10.2.1 I have three feature classes in geodatabase "MyGDB.gdb".  Let's call the feature classes "A", "B", and "C" respectively.  I am merging the feature classes using Geoprocess > merge.  Here are the exact steps:

Geoprocessing>Merge
Select feature classes "A", "B", and "C" as the input datasets
Save the output dataset as "Merged" in "MyGDB.gdb"

After creating the merged dataset, I am running the following code in a python script to count the number of data points in each feature class and also in the merged data set:

lyrfile1 = 'C:\\ArcGIS\\Python\\MyGDB.gdb\\Merged'
lyrfile2 = 'C:\\ArcGIS\\Python\\MyGDB.gdb\\A'
lyrfile3 = 'C:\\ArcGIS\\Python\\MyGDB.gdb\\B'
lyrfile4 = 'C:\\ArcGIS\\Python\\MyGDB.gdb\\C'

result1 = int(arcpy.GetCount_management(lyrfile1).getOutput(0))
result2 = int(arcpy.GetCount_management(lyrfile2).getOutput(0))
result3 = int(arcpy.GetCount_management(lyrfile3).getOutput(0))
result4 = int(arcpy.GetCount_management(lyrfile4).getOutput(0))

print result1
print result2
print result3
print result4

The results are as follows:
151983
19394
75578
57011

This is good and expected.  Since the three feature classes are merged, I would expect to see that the merged data set has the same number of data points as each individual feature class combined,, and it does (57011+75578+19394=151983).

Here is where the problem starts:  when I try to merge the datasets using arcpy.Merge_management() and then use arcpy.GetCount to count the datapoints, only a fraction of the datapoints are returned (57011).  Here is the code:

fClasses = []   #initialize array to hold feature classes

#Set counter for loop
i=0

#Loop through MyGDB.gdb and add each feature class name to array
try:
        # Get a list of the featureclasses in the gdb
    fcs = arcpy.ListFeatureClasses()

    #Iterate
    for fc in fcs:
        print fc
        fClasses.append(fc)
        i = i+1

except Exception as err:
    arcpy.AddError(err)
    print err

#Now that we have a list of all feature classes, attempt to merge them into one class
try:
    arcpy.Merge_management([fClasses], "Merged")
    print "merge:  SUCCESS."

        #Count number of points in newly created merged dataset
    result1 = int(arcpy.GetCount_management("Merged").getOutput(0))
    print result1

except Exception as err:
    arcpy.AddError(err)
    print "merge:  FAIL."
    print err

This returns only 57011 points, versus the 151983 points in the merge file creating using Geoprocessing > Merge.  What am I missing here?


texts:
I'm trying to perform a rather simple task.. but can't find a simple solution.. I'm trying to re-project a raster layer using GDAL in python. I've read this solution , and this , and this, and others. It seems that all of them assume some parameters like pixel size or xy origins but I need to get all of the parameters from the raster itself, it suppose to be as generic as possible, the only constant thing is the output coordinate system which is wgs84 (EPSG=4326).

I've also saw different tools to perform this task like gdalwrap or getransform, so I'm a bit confused.

To sum it all up, I need to re-project the input_raster from whatever coordinate system it is on to output_raster that will be in wgs84 coordinate system.

input_raster = r"C:\Users\B3.tif"

output_raster = r"C:\Users\B3_proj.tif"

Could you help me in writing a python code using gdal to perform this task?


texts: Merging convoluted raster DEMs - what is this method and why use it?

texts:
This is more a theoretical question than practical, and I'm referring to 3x3 convolution kernels with online demo and the first comment after it.

I have a raster elevation model and want to turn it into slope/gradient using convolution. The article gives me kernels to use for horizontal gradient detection and vertical gradient detection.

I use them and so have two convoluted DEMs; one with horizontal gradient and one with vertical. To get the full effect I must merge them together.

The comment underneath the article suggests that to merge these I should use the equation:

pixel = SQRT(pixel_convolved_along_x ^ 2 + pixel_convolved_along_y ^ 2)

I do that and the result looks good. In fact I compare it to a slope generated directly by a slope-creation tool and the results look identical.

So I have a process that works and a result I'm happy with - I just don't understand why that equation is used to merge cells. It looks like Pythagoras to me (a=sqrt(b^2+c^2)).

So why that and why not simply average the two cells using (a+b)/2?

Can anyone help me understand this?


texts: Drawing of building outlines (rooftops) based on low resolution aerial images

texts:
Google maps provides street and rooftop maps from corresponding aerial imagery, but not for all places.

Please is there a method or available software to help outline rooftops from a set of map-tiles, for areas not covered by Google maps?

how to convert to this?


texts: Can GpsBabel join tracks into single track?

texts:
I sometimes need to merge/join multiple tracks from multiple GPX files, and end up with a single track in a single file.

Even if the tracks are in the right order + direction, it looks like GpsBabel won't perform this task if the files have no <time> ("trackfilter-init: Found track point at [lat] without time!") or overlapping <time> ("trackfilter-pack: Tracks overlap in time!").

gpsbabel.exe -t -i gpx -f 01.gpx -i gpx -f 02.gpx -i gpx -f 03.gpx -x track,pack -o gpx -F merged.gpx

Do you confirm that it cannnot be told to simply ignore <time>, and merge all the points?

Thank you.

Edit: The solution is mentioned at the bottom of the section "Manipulate track lists (track)":

gpsbabel -i gpx -f 1.gpx -f 2.gpx -f 3.gpx -x track,merge,faketime=f20100705200000+2 -o gpx -F merged.fake.time.gpx

For some reason, it doesn't work with 19700101000000+2, ie. 1970 Jan 01, 00:00:00, increment 2s:

gpsbabel -i gpx -f 1.gpx -f 2.gpx -f 3.gpx -x track,merge,faketime=f19700101000000+2 -o gpx -F merged.fake.time.gpx
trackfilter-init: Found track point at 46.672889,-1.435020 without time!


texts: How to visualize MBTiles with OpenLayers using a PHP tile server?

texts:
This is repost as didn't get my expected answer. Please help me as I am newbie.

I have setup the PHP Tile Server 0.1 and when I am calling the PHP server it showing the following

Now I am triyng to visualize the map using openlayers. I am using the following code. I am just in the beginning stage so small mistakes may be really easy to me.

  <title>OpenLayers Example</title>

I am struggling from last night. Please help me.


texts: Identify Feature by dragging polygon over datasets using QGIS?

texts:
When using ArcGIS Desktop, I often use the Identify tool by clicking and dragging the tool over my data and identifying the features within it.

Is there a way to do this in QGIS or is setting the search radius within the map tools options the best that can be achieved for searching across multiple datasets?


texts:
I am new to GIS and still learning. It’s interesting but becoming too complicated for me. A part of my project is to determine the land cover change in some counties in Texas from 2001 to 2011.
So, a friend told me to use this approach using ArcGIS 10.3:

I downloaded the Texas NLCD data
I opened ArcMap, used the county shapefile to select the counties I’m interested in, I then overlaid the Texas land cover data (started with the 2006 data; just to try it out) on it.
After selecting the counties, I created a specific buffer zone around them.
Using ArcMap tools, I clicked spatial analyst tools →extraction→extract by mask→ it gave the option to select what land cover data I want in the input raster, I chose the Texas land cover data →then I selected the land cover data of my counties of interest into the input raster/feature mask, then selected an output file and ran  it. The result was a number of legends with different colors indicating different land cover types. I am interested in specific land cover types (four actually). What can I do to get only the legends of the land cover types I am interested in?
I also went further to click on the patch grid menu→then the spatial analyst (Fragstat interface), after which I chose the layer I wanted to use (land cover data of the counties of choice), clicked analysis by class and chose only class area under the class options menu, then I clicked run. I got an error saying Fragstat unable to run, and a table with no data output.

Is there anything I should have done? Am I using the right approach for this analysis?


texts: I need to overlay 3D models of buildings I have designed in a separate CAD software. I want to add .stl or .obj files with proper scaling on the 2D map on qgis so that I can make a 3D map with 2D base image. Is it possible in QGIS? If yes, how do I go about it? if not, what other alternatives do I have for doing this?


texts: Multiple points at same location in CartoDB - how to make clickable?

texts: I'm trying to create a simple map based on data from our website. Whenever someone wants to access a specific file, they have to fill out a form that includes their ZIP code and some other information. To display where we're getting requests from, I want to have a map with clickable dots, and a pop-up window that includes more detail about the requester. This is pretty straightforward with the "Simple points" wizard, except for when there are multiple requests from the same ZIP code. They all get geocoded to the same location and you can only click one of them. Is there any way to make this work in CartoDB?


texts: Do elevation measurements differ between Coordinate Reference Systems? If so, can they be converted?

texts:
I have global elevation data captured in relation to the WGS84 ellipsoid, and a local elevation dataset which uses the OSGB 1936 CRS. Would one expect the values to be different, given that the two reference systems use a different ellipsoid, or are elevation measurements taken in relation to some common datum independent of CRS? If the values are different, how can they be converted from one CRS to another?

This is well outside of my usual field of expertise so I apologise if this is a naive question.


texts:
I have issues extracting geometrical shaped data from a GRIB raster of soil moisture. The file is for example https://gimms.gsfc.nasa.gov/SMOS/jbolten/FAS/L03/20100113_20100115.as1.grib

Using either GDAL or Rasterio I can read the file and display it here.

Its CRS is CRS({'proj': 'stere', 'lat_0': 90, 'lat_ts': 60, 'lon_0': 1, 'k': 1, 'x_0': 0, 'y_0': 0, 'a': 6371200, 'b': 6371200, 'units': 'm', 'no_defs': True})

or with gdalinfo it reads

Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: 20141227_20141229.as1.grib
Size is 1440, 600
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371200,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",1],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (662126.628862303332426,59159218.394239947199821)
Pixel Size = (25000.000000000000000,-25000.000000000000000)
Corner Coordinates:
Upper Left  (  662126.629,59159218.394) (179d38'28.48"W, 67d16'31.51"S)
Lower Left  (  662126.629,44159218.394) (179d51'32.52"W, 59d52' 0.26"S)
Upper Right (36662126.629,59159218.394) (149d12'45.70"E, 70d36'45.13"S)
Lower Right (36662126.629,44159218.394) (141d17'58.79"E, 66d35'39.62"S)
Center      (18662126.629,51659218.394) (161d 8'15.06"E, 65d34'25.05"S)
Band 1 Block=1440x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC (Ground or water surface)
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=undefined [-]
    GRIB_ELEMENT=var255
    GRIB_FORECAST_SECONDS=86400 sec
    GRIB_REF_TIME=  1419638400 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[-]
    GRIB_VALID_TIME=  1419724800 sec UTC

The issue occurs if I want to extract a region from the raster, using a predefined shape.

Using rasterio.warp.transform_geom with source CRS EPSG:4326 and dest CRS of the GRIB file above, I end up having a huge shape and of course the new one does not intersect with the raster.

Could anyone show me how to do it right with this GRIB format? I have done such operation without any problem on netcdf, hdf, and geotiff files.


texts: Difference between S3GeoTiffReader and GeoTiffReader?

texts:
What's the difference between S3GeoTiffReader and GeoTiffReader provided by Geotools library?
I need to process a GeoTiff image retrieved from AWS S3 bucket.


texts:
I need to generate isochrones for emergency vehicles. Is there any possibility of QGIS to do this?
I checked TravelTime and QNEAT plugins but unfortunately, I could not find any way.
With TravelTime there is not any travel mode considering emergency vehicle speed in TravelTime plugin. Since emergency vehicles don't travel in ordinary speed limitations I need to configure the travel speed.
Do you know any other way, not only using QGIS, for generating this type of isochrone?


texts: Does raster snapping involve a transformation?

texts:
Arc Desktop 10.5.1

Background

I am converting an Arc Grid raster to a tif using the Raster To Other Format (multiple) tool.

In the tool's Environments... settings, I am snapping the output tif to a different tif raster.

All rasters involved have identical cell sizes (three feet), and all have the same horizontal and vertical CRS/projections.

My question:

The Environments... help implies that only an x/y shift occurs with raster snapping, without any sort of transformation (such as in a reprojection).

However, the Environments... help also states, within the Raster Storage section:

“Resampling is the process of interpolating the pixel values while transforming your raster dataset.  This is used when the input and output do not line up exactly, when the pixel size changes, when the data is shifted, or a combination of these.”  My emphasis.

My understanding has long been that a snap by itself would not invoke a transform (and thus resampling).  But after a closer look at the above help (IMHO ambiguous) I'm not so sure!  I need to know if it does so that I can then decide what resampling technique I want to use, so I'm asking for assistance: does a raster snap involve transforming, and thus resampling?


texts:
Today I've been messing around with QGIS installations and I think I've accidentally caused a problem: QGIS is no longer able to execute the GDAL raster tool 'Fill no data' and other GDAL tools that use Python.

I am using Windows 10 and installed QGIS 3.10.1 from www.qgis.org using the standalone QGIS version. Later, I uninstalled that installation from my pc, and installed QGIS and GDAL and GRASS GIS using the OSGeo4W Network installer (also through www.qgis.org). Then I couldn't figure out how to remove the OSGeo4W installation, so I just deleted the folder 'OSGeo4W' on my C:\ drive. And then I installed the standalone QGIS version again.

At this point, whenever I try to execute the GDAL tool 'Fill no data' on a raster layer in QGIS, I get the following error. Other GDAL tools that use Python generate the same error.

2019-12-11T17:20:16     INFO    python3 -m gdal_fillnodata -md 10 -b 1 -of GTiff C:/Users/username/Downloads/ahn_ruw/M_32GZ2.TIF C:/Users/username/Downloads/ahn_fill/test.tif
2019-12-11T17:20:16     INFO    GDAL execution console output
             C:\PROGRA~1\QGIS3~1.10\bin\python3.exe: No module named gdal_fillnodata

What I think is happening is that my current QGIS installation is looking for a non-existent Python/GDAL installation in the old OSGeo4W folder. Or some other confusing with different installations of Python/GDAL. I think it is a Python problem.

I have tried this but that didn't work. I've also tried changing the system variable PATH in QGIS under Settings > System to the folder C:\Program Files\QGIS 3.10\apps\Python37 but that didn't work. It seems that the correct path (C:\Program Files\QGIS 3.10\apps\Python37) is already in the system variable PATH.

Does somebody know how to fix this?


texts:
I have discovered that cloudmade tiles service has been discontinued, so to switch to mapbox tiles, i imagine that i can find a gallery (like in cloudmade ) where i can choose the tiles style convenient. i searched a lot, i've just found a demo made by a techslides.

where i could find a place where users shares thier style or any gallery?


texts: What are the parameters for my connection QGIS

texts:
I am new QGIS, I was wondering what the parameters are for the New MSSQL Connection in the dialog window?

Provider/DNS: ? when local ? tried leaving empty

Host:         ? when local ? tried my computer user name/server instance name

Database:     full path to file I assume

Error:


texts: How to remove complicated shapes in ArcGIS - generalization

texts:
What can be the way of simplifying polygons in ArcMap? I have a shapefile with buildings and I have to make them as simple as it's possible. I used already some tools, like: dissolve, union, generalize, integrate and even buffer. I still have something like this to remove:

Small angle between buildings makes problems with mesh generation I want to do afterwards. Any ideas how to remove it massively for all objects? I know I can increase buffer distance, but maybe there is something which will not change geometry that much.

EDIT:
After using Simplify Polygon or Simplify Bulilding some areas are simplified (left side of picture) but not those I would be interested in:

And as an answer for the duplicate - my question is different, I don't want to delete overlapped lines, just remove holes which have very small angles.


texts:
I have made some processing on some satellite images as numpy arrays. I want to export the data as tif images again to see the result I read ùy image with rasterio
sentinel2 = rasterio.open('sentinel2.tif')

and my processed array is date1seg
According to this >>> new_dataset.write(Z, 1) I used that to export the image but this error appears
sentinel2.write(date1seg, 1)
AttributeError: 'DatasetReader' object has no attribute 'write'


texts: What's the best way to serve a changing set of Geotiffs as a WMS layer?

texts:
I have a set of geotiffs that will be continually updating and I'm trying to find the best way of automating this process. At any point in time, I want to be able to make a WMS request against a single layer name that returns the composite single jpeg etc. and I'd also like to be to constrain that WMS search using time fields. I am using the OpenGEO suite with Geoserver and a PostGres DB.

I envisaged having a folder on my server system that I would transfer each new geotiff as it comes in and probably running a curl command to somehow add that to the existing collection. Looking around, it was unclear what would be the best approach that would support this.

It seemed like I can use the Geoserver REST interface to add a Geotiff as a layer but it wasn't obvious if I could somehow keep on adding each new geotiff to this layer or if it was only for a 1:1 mapping of geotiff to layer name (but then maybe using a layer group and adding each layer to this group)?
Coverage stores - (Programming Geoserver 2.0.2 to add a raster data store and layer without the UI) Underdark's comment sounded promising but it was not obvious what the constraints are with using a coverage store (ongoing reading around on this one currently ;-) )
Somehow using the PostGRES to store the metadata etc. (and maybe a raster quicklook) but this approach seems more relevant for vector info rather than raster.

Currently my geotiffs are only about 35MB each but I'd like to extend this system to support multiMeg/Gig satellite imagery so I'd need to be able to support some form of pyramiding in the future so thoughts on expandability of the above challenge to support this functionality would be useful to hear (e.g. better to create multiple individual files for each pyramid level and assess which filename to access (expected to be more efficient for memory management??), OR creating a single composited pyramid file with 1x1, 2x2, etc. levels (or maybe 1x1 stays separate and the pyramid has 2x2, 4x4 etc. only))


texts: Batch Assign Coded Domain Values with Multiple SubTypes?

texts:
I have a File GeoDatabase (v10.0) with multiple layers.  Each layer has a range of Subtypes (1 to 20+).  I have created coded domain values that I would like to assign to their appropriate field.  I know the manual approach, but I am hoping to find a more automated system.  Does anyone know of a way, tool, plugin, etc. to batch assign the coded domain values to fields that have subtypes?

Attached is a graphic of the properties of one of my layers.  In the Subtypes section you can see that there are many Codes (205, 206, 207, 224, etc.).  In the Default Values and Domains section there are four (4) fields of interest (TypeComp, Type2Comp, Layer, and UserName).  In order to have the correct Domains associated to the correct field, I need to associate the Domain to the field for each Subtype code.  I want to know if there is an automated way to do this, as I have 13 layers, and each layer has multiple Subtype codes.


texts: Get distance from point to nearest edge of polygon with PyQGIS

texts:
I have 2 .shp layers in my workspace, polyLayer and pointLayer

Using PyQGIS I'd like to be able to get the distance from the point to the nearest boundary of the polygon.

I've done this sort of thing successfully in PostGIS. It seems like it must be a very commonly used query, but I can't find anything on performing such queries on .shp files using PyQGIS.

Does anyone know how this can be achieved?


texts:
I have two sets of data:

list of parcel numbers and locations which I am interested in
Address database for the county (shapefile)

In dataset 2 there can be multiple instances of the same parcel number with different addresses.

However in dataset 1 there is only 1 unique instance of each parcel number.

I would like to use dataset 1 to extract ALL instances of the same APN in dataset 2.

I have tried a table join but this will not work since table join looks for one and only one corresponding value for parcel number.

I need all corresponding values.


texts: Rasterize vector while adjusting extent and resolution to reference raster?

texts:
I have a vector file which I want to rasterize and adjust its resolution to the reference map, in QGIS.

I tried spatial_sync_raster function in R to set all the rasters into the same resolution as that of the reference raster layer. It worked but the result was not as good as expected, because the outputs came out with different resolutions (see examples below).

I want to set the resolution of "u" as that of "g".

u
class       : RasterLayer
dimensions  : 130, 228, 29640  (nrow, ncol, ncell)
resolution  : 0.01213899, 0.01213899  (x, y)
extent      : 12.111, 14.87869, 51.68928, 53.26734  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

g
class       : RasterLayer
dimensions  : 189, 332, 62748  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 12.10837, 14.87503, 51.7, 53.275  (xmin, xmax, ymin, ymax)

Does anyone know how to adjust the resolution whilst rasterizing the vector file?


texts:
I've got back into Python after a long break. I'm creating a simple data validation script at the moment. I am just trying to return the field type of a field input by the user. I'm using arcpy.AddMessage to validate my steps as I go along. I've tried googling and looking at multiple questions on here and have found nothing that works. Here is my code:

    #VARIABLES
gulyShp = arcpy.GetParameterAsText (0)
#GULLY SHAPEFILE
sDep = arcpy.GetParameterAsText (1)

arcpy.AddMessage(sDep)

sDTy = arcpy.Field(sDep.Type)
arcpy.AddMessage(sDTy)

All I want is the type of the field sDep is assigned to be stored in a variable.


texts: Result view of a multiple polygons query doesn't fit to all of them and is centred just in one

texts:
I have built a Mapbox web app with some polygon with different attributes.
The problem is when I query a table and the result is that many polygons have the same attributes. The zoom  centre in one of the polygons, but doesn't show all the polygons in the zoom view, with this common attribute as I expected.
So at some point the query below doesn't provide the correct information to create the bounding boxes I guess, but I cannot manage to see what is failing if that is the error.
This is the query I used to retrieve the polygons depending on the  title_no with other additional data from additional tables:
SELECT
            pt.title_no,
            num_polygons::integer,
            rec_status,
            round(st_area(pt.geom)) area,
            array_to_json(ARRAY[
            st_x(st_transform(st_pointonsurface(pt.geom),4326)),
            st_y(st_transform(st_pointonsurface(pt.geom),4326))
            ]) center,
            ARRAY( select json_array_elements_text((st_asgeojson(st_envelope(st_transform(pt.geom, 4326)),15,1)::json->>'bbox')::json))::float[] bbox,
            st_asgeojson(st_transform(pt.geom, 4326)) geom
        FROM polygon_table pt
      ;

I don't know if it is related but also happens with large polygons doesn't show the whole area and needs to be zoom to see it.
EXTENDED QUESTION:
That is how I add the fitBounds to the GIS app,
map = new mapboxgl.Map({
      container: container.current,
      style: 'url to map box style'
      center: [-2, 53],
      zoom: 5.75,
      accessToken: mapAccessToken,
      antialias: true
    });

    setMap(map);

    map.on('load', () => map.fitBounds(boundingBox, { duration: 0, maxZoom: 17, padding: 30 }));


texts:
I am wondering why ogr2ogr is dropping certain features when exporting to GeoJSON and how to change this behavior.

In the following picture, all the features extracted by ogr2ogr are coded in fuchsia. Almost all of the forest is missing.

But when I load the raw OSM extract in JOSM, we see the forest is represented, though its relation is truncated at the northern border of the bounding box.

I have verified that the key-value pair in the OSM xml is landuse -> forest, and I have configured my customized osmconf.ini to include the landuse key for the multipolygon attributes.

So why does ogr2ogr exclude the truncated forest relation from the multipolygon GeoJSON? And is there a way to get it to include the truncated relations?

I should note that there are other instances where ogr2ogr does generate multipoloygons where the original feature has been split along a bounding line.

EDIT:
Here is a link to the region:
https://www.openstreetmap.org/#map=12/39.0835/-123.2851

And this is the relation in question:
https://www.openstreetmap.org/relation/1015303


texts:
I'm writing a QGIS plugin and I need to create a raster layer as the average of other rasters. How can I do it? I have tried to use gdal_calc.py through os.system, something like:

cmd = 'gdal_calc.py
       -A first_layer.tif
       -B second_layer.tif
       --outfile=AVG.tif
       --cal="(A+B)/2"'
os.system(cmd)

but it seems that the command (which works perfectly fine in the command line) is not executed from my QGIS plugin. What should I use to create my average raster?

I finally found a solution, and I share it for those who might be interested:

from gdal_calculations import Dataset, Env

def aggregate_layers(self, layers):
    datasets = []
    for l in layers:
        datasets.append(Dataset(l))
    sum = datasets[0]
    for i in range(1, len(datasets) - 1):
        sum += datasets[i]
    Env.overwrite = True
    avg = sum
    avg /= len(datasets)
    root = '/home/kalimaha'
    file_name = root + '/average.tif'
    avg.save(file_name)
    if self.dlg.open_in_qgis.isChecked() is True:
        rl = self.iface.addRasterLayer(file_name, 'Average Layer')

I have tried @xunilk's code, but my QGIS simply crashes. I had to edit it a bit, this is the final code:

def aggregate_layers(self, layers):

    entries = []
    for l in layers:
        fileInfo = QFileInfo(l)
        baseName = fileInfo.baseName()
        rlayer = QgsRasterLayer(l, baseName)
        QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')
        layer = QgsRasterCalculatorEntry()
        layer.ref = rlayer.name() +'@1'
        layer.raster = rlayer
        layer.bandNumber = 1
        entries.append(layer)

    for e in entries:
        QgsMessageLog.logMessage('e: ' + str(e), 'Geobricks TRMM')

    fileInfo = QFileInfo(layers[0])
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(layers[0], baseName)
    QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')

    expression = '(' + entries[1].ref + ' + ' + entries[0].ref + ') / 2'
    QgsMessageLog.logMessage('Expression: ' + str(expression), 'Geobricks TRMM')
    calc = QgsRasterCalculator(expression,
                               '/Development/GIS/average.tif',
                               'GTiff',
                               rlayer.extent(),
                               rlayer.width(),
                               rlayer.height(),
                               entries)
    QgsMessageLog.logMessage('calc: ' + str(calc), 'Geobricks TRMM')

My layers are not in QGIS yet, so I have to load it from the filesystem, the array nemd layers contains the absolute paths of the layers I need to aggregate. So, when I do calc = QgsRasterCalculator(... QGIS crashes.


texts:
I am using MOD09Q1 now and I would like to calculating NDVI whole year data. One year, I have 46 band 1 images and 46 band 2 images.

list.files

I creating list of band 1 and list of band 2.

Band1list  <- list.files(path = ".", pattern = "*.b01.tif", all.files = FALSE,
                         full.names = FALSE, recursive = FALSE,
                         ignore.case = FALSE)

Band2list  <- list.files(path = ".", pattern = "*.b02.tif", all.files = FALSE,
                         full.names = FALSE, recursive = FALSE,
                         ignore.case = FALSE)

raster::stack

Stacking all 46 band 1 (band 2) images to 'StackBand1' ('StackBand2')

StackBand1 <- stack(Band1list)

StackBand2 <- stack(Band2list)

NDVI calculation via for loop

I try to stack images of same date (e.g. "2012.01.01.Smb01.tif" and "2012.01.01.Smb02.tif") and calculate NDVI

for (i in 1:length(Band1list))
{
  StackBand1[[i]] == i

  StackBand2[[i]] == i

  Stack <- brick(StackBand1[[i]],StackBand2[[i]])
}

In here, I could not execute loop for all 46 images. It can be run all and result is showed for the last couple only.


texts: Transforming polygon from Arc 1960 UTM Zone 37S to WGS 1984 UTM Zone 37S?

texts:
I have a polygon which is in Arc 1960 UTM Zone 37S coordinate system and I need to transform it to to WGS 1984 UTM Zone 37S. I have not been successful in doing this using ArcMap.

Is there a solution for this?


texts: Order/Sorting of Iterate Field Values (ModelBuilder)

texts:
I am using Iterate Field Values to Iterate my Model through a list of Strings (DBF-File). However, the Field Values-Iterator will always sort my Field by alphabet (ascending). It does not matter how the rows in the *.dbf are ordered.

Is it possible to run the model using a different sorting method?

Example of dbf-file:

house
mill
water
tree
cafe
garden

.. will run as

cafe
garden
house
mill
tree
water

I would like to run the Model by the first sequence of strings, not the second one. Is this possible?
I tried adding a second column with numbers from 1 up to x, but the Iterator ignores this field.


texts: Point-Line shortest distance and corresponding angle using FME

texts: I want to determine the shortest distance and especially the corresponding angle from a point to a polyline. A solution using FME would be strongly preferred.


texts: Leaflet markers at same position: dynamically display all markers infos through one marker only

texts:
I created a map using Leaflet, with a filtering system, custom icons with the leaflet.awesome-markers plugin, and jquery functions so that when you click on the markers, the corresponding information, inherited from a geoJson file, appears in turn in a div. (Every time I click on a marker, the marker info appears in my div, dynamically). I also use the markercluster plugin to display clusters when several points are in the same place.

Here is the rendering:

The filters are checkbox (multiple) so user can click on all filters to see all markers.

The markers represent hospitals, so when there is a cluster, it shows multiple markers for the same hospital.

The clusters are working well, and if I click on it, I can see the clustured markers.
Here I have two markers for the same hospital. And that's my problem!

Now I have to change the behavior I've gotten.

What I would need  to do:

Keep the cluster system, but make sure that there is only one marker when zooming, and that the information of all the markers of the cluster are displayed in my dynamic div. (with a tabs system or other?)
So I want to keep the cluster system to show that there are multiple informations at that point, but show only one marker when clicking on the cluster.

Is it possible ?

EDIT 1:
I edit this thread because after many search, I read https://digital-geography.com/working-with-clusters-in-leaflet-increasing-useability/ and this one and saw it was possible !

So Now just after my getJson function I added

var promise = $.getJSON("examen.json");
var clusters = L.markerClusterGroup( { <-- this is new
    spiderfyOnMaxZoom: false,
    showCoverageOnHover: false,
    zoomToBoundsOnClick: true
} );

Then, at the end of my onEachFeature function, I added

clusters.on('clusterclick', function(a){
    var markers = a.layer.getAllChildMarkers();
    for (var i = 0; i < markers.length; i++) {
   $('#layer_infos .fill').append(html);
    }
});

It is working partialy.

In the documentation, it is written

  Returns an array of all markers contained within this cluster

But for me it doesn't show "this" marker info, but all the markers of my geoJson file... ! (image below with the scroll bar)  maybe because my cluster function is wrapping all my markers in the code ? So how can I limit it to show only the clicked marker info ?


texts:
I am looking for a non-manual way to stretch a colour range over two raster layers.
Lets say I have two raster with data showing mean annual precipitation. In order to compare them visually I would like to stretch a colour range over both layers.
So in my case it should range from 320 - 900 and for example range from light blue to dark blue over both layers.


texts: Routing two not-node points and road side with pgrouting

texts:
I'm starting to use pgrouting and I have a couple of problems. In my project I have a lot of points which are near one of the road sides.
I'd like to know how I can route two of these points since they aren't nodes of the topology. This route have to consider the road side in which they are and finish in this side and not in the other one.


texts: Projection of set of point on a line with chainage

texts:
I am new to QGIS and in fact new to GIS mapping as whole.

I have a line in shp format (which does not have chainage & I want to add chainage added to it).

Chainage is the distance from the starting point of a line. I have tried using Qchainage Plugin to extract the distance (chainage) from one of the starting point. After that I used NNJoin Plugin to join two layer. Problem is, Qchainage doesn't add the attribute to the same layer of line which could be related to distance (chainage) & then related to multipoints.

Then I have several datapoints in the vicinity of the line. Attribute tables have been provided to all datapoints except the chainage along the line.

It would be helpful to first add chainage to the line & then project these points on the line with an added attribute termed as chainage.


texts: Retrieve tile path of superoverlay in Google Earth Plugin using Javascript

texts:
I've created a superoverlay using MapTiler and am displaying that in a webpage with the Google Earth plugin.  I'd like to identify the current tile (level, x, y) being displayed with a mousemove event.  I'm able to get the current lat/lon using getGlobe(), e.g.

google.earth.addEventListener(ge.getGlobe(), 'mousemove', function(e) {
    var lat = e.getLatitude();
    var lon = e.getLongitude();
}

The biggest hurdle seems to be getting the correct superoverlay level based on camera altitude as the transition between levels is dictated by the minLodPixels parameter in the KML files.  I'm not sure how best to determine when a tile is, say, 128 pixels wide on the screen, especially when this will be based on how large the map window div is.

Once the level is determined, it should be trivial to figure out the X, Y offsets of the tile.  I'm also assuming that the minLodPixels/maxLodPixels parameters will stay consistent across all tiles in the superoverlay (MapTiler has these at 128/-1).

Any thoughts?

Update 2012.02.22

I've delved a bit further into this and found some good documentation on minLod in the KML documentation.  From http://code.google.com/apis/kml/documentation/kmlreference.html#region,

  Google Earth calculates the size of the Region when projected onto screen space. Then it computes the square root of the Region's area (if, for example, the Region is square and the viewpoint is directly above the Region, and the Region is not tilted, this measurement is equal to the width of the projected Region). If this measurement falls within the limits defined by  and  (and if the  is in view), the Region is active.

So if I'm able to find out the screen real estate of each hierarchical tile at a particular LookAt, I should be able to find the level I'm at based on what tile minimally exceeds the minLod threshold.  As a test, I've read the SuperOverlay's root KML ('doc.kml') and retrieved the coarsest tile bounds:

var INFO = '/path/to/doc.kml';
$.get(INFO, function(data) {
    // Get the corners of the coarsest resolution
    var box = $(data).find("LatLonAltBox");
    north = parseFloat($(box).find("north").text());
    south = parseFloat($(box).find("south").text());
    east = parseFloat($(box).find("east").text());
    west = parseFloat($(box).find("west").text());
});

From there, I calculate area on the mousemove event of GEWindow.

google.earth.addEventListener(ge.getWindow(), 'mousemove', function(e) {

    // Track the screen coordinates of the coarsest tile
    var view = ge.getView();
    var alt_mode = ge.ALTITUDE_RELATIVE_TO_GROUND;
    var ne = view.project(north, east, 0.0, alt_mode);
    var se = view.project(south, east, 0.0, alt_mode);
    var sw = view.project(south, west, 0.0, alt_mode);
    var nw = view.project(north, west, 0.0, alt_mode);

    // Get the square-root of area (using equation 3 from
    // http://mathworld.wolfram.com/Quadrilateral.html)
    var qX = (se.getX() - nw.getX());
    var qY = (se.getY() - nw.getY());
    var pX = (sw.getX() - ne.getX());
    var pY = (sw.getY() - ne.getY());
    var area = Math.sqrt(0.5 * Math.abs(qX * pY - pX * qY));
});

When I track the tile's visibility in GE though, it disappears at a value higher than 128 (the minLod) when viewed without tilt, and appears at a value lower than 128 when viewed with tilt.  I don't have enough rep to post pictures, but I've posted example screenshots here.  In these pictures, the white polyline shows the tile boundary.  I'm obviously still doing something incorrectly.


texts: Is it possible to redesign OpenStreetMap tiles?

texts:
We have an application working with OpenLayers and are using OpenStreetmap for the map tiles, however they don't look great.

Is there anyway to style these maps?

I read something about Mapbox but it uses Leaflet API, however we use OpenLayers and our app is heavily using it.  What options do I have?


texts: Contouring range only using ArcGIS Spatial Analyst?

texts:
I have a grid (in .img format) with  pixel values ranging from -10,000 to 10,000, which i would like to contour. but i would like to contour only -3000 to 0 range. is this possible to set this up during contouring? or is there a way i can restrict the raster grid to display only -3000 to 0 before i start contouring?

Using Arc 9.3.1, with one extension (spatial analyst), nothing else.


texts:
When zooming into my point layer it's just one big blob of dots. I want to be able to zoom into see individual points but scale 1:2 goes nowhere near to this extent.

The layer was generated in ArcMap in GCS_WGS_1984. When I published the layer I used EPSG:32610

I computed the native bounding box from the data and Lat/Lon box from native bounds.


texts: Stop Virtual layer from automatically recalculation

texts:
QGIS virtual layer has big performance impact. Even when panning or zooming in the canvas, a virtual layer could take seconds to minutes to refresh. From my speculation, the reason might be the query it has to go through.
My questions is, would that be possible to cache a virtual layer, so it doesn't have to go through all the query until a user asks for it?
Say, I may update the database once a day, so I would like the virtual layer to update also once a day on my command, instead of refreshing every time I zoom in or change a layer visibility.
Specifically, I have a virtual layer lot_comparison combining two layers, existing_lot(containing 900 records) and planning_lot (containing 300 records).
the query goes like this,
select st_union(pl.geometry), el.*, pl.*
from "existing_lot" as el, "planning_lot" as pl
where el.fid = pl.fid
and (el."lotIDtype" in ('use previous lot index' ,'need subdivision')
or pl."lotIDtype" = 'need new lot ID')
group by pl."identifier"


texts:
I am trying to install pyproj. I downloaded the zip folder on this page pyproj.

Then I copied it to C:\Python27\ArcGIS10.3\Lib. However, when I try to import the module, I get the message:

  ImportError: No module named pyproj.

I also tried to run the setup.py, but I get the message ImportError: No module named setuptools. I read somewher that it is working when I install Basemap, which is part of Matplotlib, but I cannot find pyproj there.

Any help?


texts: Managing Committed RAM % in ArcGIS Server

texts:
I have a server with 16GB of RAM, and it's holding roughly 156 Image Services.

Managing the % of committed bytes in use, I am using roughly 83%. Or about 13.28 GB.

The server is now kicking out some of the Image Services, and when I try to restart/add them back, I get errors.

Does the server really need 2.7 GB of breathing room?

Are there ways to get around this without using 2 CPUs or increasing RAM again?


texts: Leaflet - Opacity method for L.tileLayer.wms?

texts:
I want to use opacity on one of my WMS layers in Leaflet. The documentation has an opacity option for L.tileLayer but not explicitly for L.tileLayer.wms, I only see a transparent option (not supported for the WMS I want to use, I'm using a cascading WMS and it defaults to gif I think, I have not found a way how I can cascade it as a PNG). Does this mean that you can't use opacity only on L.tileLayer and not for L.tileLayer.wms?

This does not do anything:

var gemini = L.tileLayer.wms('rootUrl', {
    layers: 'myLayer',
    format: 'image/png',
    maxZoom: 21,
    opacity: 0.5;
});


texts: How to get only Intersected polygon in Open Layers 3

texts:
I am trying to draw a line on polygon which I want to split,
Example:
I want to split 616 polygon as show in image

with below code when I tied to get intersected polygon I am getting 616,617,615 and Roads adjacent to 616.But I want only 616 which cut by line

draw.on('drawend', function (evt) {
                    var plcd = [], featureGeom;
                    var justNowFeature = evt.feature;
                    var line = evt.feature.getGeometry();
                     featureGeom = justNowFeature.getGeometry().getCoordinates();

                    console.log("FEATURESGEOJSON  " + featureGeom);
                    var features = plotboundary_layer.getSource().getFeatures();
                    for (var i = 0; i < features.length; i++) {
                        if (line.intersectsExtent(features[i].getGeometry().getExtent())) {
                            plcd.push(features[i].values_.plcd);
                        }
                    }
                });


texts: As many probably know, the MacOS installer for QGIS is not as turnkey as the windows version.  You have to basically install 3-4 sets of packages, only one of which end up in the Mac apps directory (e.g., GDAL) . If I want to completely uninstall QGIS, how do I go about uninstalling the packages like GDAL?


texts: PostGIS geometry intersection test results different for Point and Linestring

texts:
Can anyone please tell me what I'm doing wrong, and how do I avoid the problem? I'll attempt to clarify what I'm doing as things go along.

The geom column is a multipolygon column, all values are SRID=4326.

archive=> select st_intersects( GeomFromewkt('SRID=4326;LINESTRING(-16.899216666666668 12.610816666666667,-16.899216666666668 12.610816666666667)')
, z.geom) as intersects
from zone z
where z.zone_id=2;
 intersects
------------
 f
(1 row)

archive=> select   st_intersects( GeomFromewkt('SRID=4326;POINT(-16.899216666666668 12.610816666666667)')
, z.geom) as intersects
from zone z
where z.zone_id=2;
 intersects
------------
 t
(1 row)

This is the text representation of the geom:

MULTIPOLYGON(((-18.7273147400488 12.8862644964914,-18.7334078119781 12.3609663554343,-16.5597662148828 12.3796303212717,-16.6302631480867 13.1209501528692,-18.7273147400488 12.8862644964914)))

How can the linestring not intersect with the geometry, yet a single point does? Is this some sort of rounding error, or a problem because the linestring has a length of zero? If the length is the problem, how do I avoid this problem?


texts: Is there a way to change ST_asGeoJSON return values?

texts:
The function ST_asGeoJSON returns something like:

"st_asgeojson":
   {
     "type":"MultiLineString",
     "coordinates":[
        [
          [-3.13512629299993,39.7174353130001],
          [-19.13510345099996,2.7174767110001]
        ]
      ]
    }

And I wonder if there's a way to simplify or customize the JSON it returns, as I only need the coordinates, to something like:

{
  coordinates: [
    [-9.13512629299993,38.7174353130001],
    [-9.13510345099996,38.7174767110001]
  ]
}


texts: How to represent microbuses (hail and ride transport) in GTFS (for OpenTripPlanner server)

texts: I am doing transit database for my city in GTFS format. And I have an problem with microbuses (hail-and-ride transport), because this transport doesn't have stops. As a workaround, I can made stops every 100-200 meters, but maybe it's possible to do easier.


texts: Split By Attributes creates no Shapefiles?

texts:
I want to create for each polygon of a shapefile a separate shapefile with the Split By Attributes tool.
Input is a polygon shapefile and as split field I choose the ID.
This tool works fine but in the results I cant find any shapefile, only .dbf , .dbf.xml and .cpg.
But why is there no .shp?

As output I tried a workspace and a geodatabase, both didn't worked.


texts:
I'm new to QGIS and work on an air-gapped system where I have access only to QGIS 2.14.9-Essen. I create a WMS layer, which "works" but performance is abysmal. So I enable caching (Settings > Network) with a large disk space.

If I simply pan back and forth by a very small amount (1-2 pixels), I expect it to use the cached files. Instead it apparently fetches new files on each pan, as indicated by the progress bar and a half-second delay. The number of files in the cache area continues to grow with each micro-pan. QGIS is caching ever-more tiles but never using them.

It might help if I understood the file structure. Under the cache area, QGIS created sub-directories data1-data8, and the cached tiles are distributed among those sub-directories, with filenames like .../cache/data6/i8duj3mh.d. What is a *.d file and why the sub-directories?


texts: Writing a blob with Esri.FileGDB.Row.SetBinary return an empty value when reread the blob field

texts:
I am trying to write an image/png in an ESRI FileGDB field.

I use the Esri.FileGDB.API 1.4, both 32/64, in C#.

The type of the field is "esriFieldTypeBlob".

I am currently able to write the image, without errors, as follows:

Table table = gdb.OpenTable(tablename);
Row newRow = table.CreateRowObject();
byte[] fba = null;
System.Drawing.Image img = System.Drawing.Image.FromFile(imgpath);
using (MemoryStream ms = new MemoryStream()) {
    img.Save(ms, System.Drawing.Imaging.ImageFormat.Png);
    fba = ms.ToArray();
}
ByteArray ba = new ByteArray((uint) fba.Length);
ba.byteArray = (byte[]) fba.Clone();
newRow.SetBinary(fieldname, ba);
table.Insert(newRow);

But when I re-read the field using GetBinary, the ByteArray is 0 length.

What am I doing wrong?

Has anyone had the same problem?


texts: Getting pictures as popup images on spatial points in QGIS?

texts:
The attached images below show a road route (symbolised as a grey dashed line) and point locations of graffiti markings on public properties - these points are symbolised as a black cross-hair.

I have taken pictures of these graffiti markings and linked their path locations of where the .jpg pictures are stored to the point coordinates (see image below).

It is unclear as to how I can get these images to appear as a popup window upon clicking these points as I'm having difficultly following the steps highlighted in previous posts:

Show images related to features in QGIS?
How to get image pop ups in QGIS?

Is there a simple way of doing this in the latest version of QGIS (3.2.3)?


texts: How do I download tracks from the GARMIN GPSMAP 78 usin DNRGarmin?

texts:
I have been able to download waypoints to DNRGarmin by accessing the GPX file. However, I have only been able to view waypoints.

How can I download tracks?


texts: Openlayers custom control div position and visibility does not appear on map

texts:
I want to createa a sample custom button to appear on my openlayers map.

     <button class="ol-zoom-in" type="button" title="Show Me">O</button>

And my css is like this

.map {
    width: 100%;
    height: 300px;
}
.ol-custom{
  top: 4.5em;
  left: .5em;
}

But my custom button does not appears. I used openlayers css classes that used for other controls but did not work.

Demo


texts: How can I change mouse icon on hover a WMS layer feature using OpenLayers and GeoServer

texts:
My unfinished app is:

"http://190.5.81.197/MAPA_OL/mapa/matricula.html#"

I just want to know why I can't change mouse icon when the mouse is over a marker on a WMS layer. I am using GeoServer to create the WMS layer.


texts:
I have a huge interrogation regarding the performance difference between Rtree  and STRtree included in Shapely. Below is a short code where we create randomly 250 000 triangles (defined by 4 vertice) in a plane of 10000 by 10000.  We also create 10000 square lines (5 vertice) and for each one, we try to find the triangles each one intersects. we test both Rtree and STRtree idependently.  What surprise me is that Rtree is almost 100 times faster to find the intersection (not counting of course the time to load the features into RTree but you do it only one- around 20 seconds).  See the result below.  One good point both trees are finding the same number of potential intersections.

My question: Is it normal that RTree is more than 100 times faster than STRtree? Did I code something wrong somewhere? Is it the way STRtree is implemented with Python overhead? Something else?

Here's the code:

from shapely.geometry import LineString
from shapely.strtree import STRtree
import random, time
from rtree import index

lst_lines = []
lst_intersects = []

# Create the triangles
for i in range(25000):
    x = random.random() * 10000.
    y = random.random() * 10000.
    coords = [(x,y),(x+5, y+5),(x,y+10),(x,y)]
    lst_lines.append(LineString(coords))

# Create the bounding boxes
for i in range(10000):
    x = random.random() * 10000.
    y = random.random() * 10000.
    coords = [(x,y),(x+15,y),(x+15,y+15),(x,y+15),(x,y)]
    lst_intersects.append(LineString(coords))

# Create shapely STRtree
tree = STRtree(lst_lines)

# Create RTree
idx = index.Index()
for i, line in enumerate(lst_lines):
    idx.insert(i, line.bounds)
print (time.time())

sec1 = time.time()

# Find the intersection with STRtree
str_tree_nbr = 0
for intersect in lst_intersects:
    str_tree = tree.query(intersect)
    str_tree_nbr += len(str_tree)

sec2 = time.time()
print("Seconds since epoch =", sec2-sec1)
print ("STRtree number: ", str_tree_nbr)

# Find the intersections with RTree
rtree_nbr = 0
for intersect in lst_intersects:
    rtree = idx.intersection(intersect.bounds)
    rtree_nbr += len(list(rtree))

sec3 = time.time()
print("Seconds since epoch =", sec3-sec2)
print ("Rtree number: ", rtree_nbr)

Results:

    Seconds for STRtree = 120.86289858818054
    STRtree number:  12622
    Seconds for RTree = 0.8911409378051758
    Rtree number:  12622

The result for the intersection time is also not linear it grows  exponentially when I add more triangles.

System used: Windows 10, conda: 4.7.12, Shapely: 1.7.0 and Rtree: 0.9.4


texts: QGIS Plugin - Automatically update field names combo box using data from another combo box

texts:
I am creating my very first QGIS plugin using two combo boxes - one to select a layer and the other to select a field.

Whilst the layer selection feature seems to work, the field name combo box is not updated when a new layer is selected.

My code is as follows...

    self.dlg.comboBox.clear()
    self.dlg.comboBox_2.clear()

    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
    self.dlg.comboBox.addItems(layer_list)
    selectedLayerIndex = self.dlg.comboBox.currentIndex()
    selectedLayer = layers[selectedLayerIndex]

    fields = [field.name() for field in selectedLayer.pendingFields() ]
    self.dlg.comboBox_2.addItems(fields)

Does anyone happen to know how I can accomplish this?


texts: Converting geotiff to kml with color gradient, GDAL

texts:
I am trying to convert a geoTIFF to a kml with gdal_translate -of KMLSUPEROVERLAY, however I don't know how to do so with a color gradient.

My geotiff is an 8bit raster with valid values ranging from 0-100. I'd like a simple output where 0-10 is gray, and 11-100 has a one-color gradient from light to dark (say, purple). This could be in classes of 11-20, 21-30, etc. or whatever, as long as it goes from light to dark.

Is this possible?

Is there a way to supply a color table with values and their corresponding RGB values? I see some options on gdal_translate but don't know where to start.


texts: Converting overlapping multipart lines into polygons in ArcGIS Pro

texts:
I have a number of multipart lines that form rectangles, each with a unique identifier. Imagine having building footprints as lines, but many intersecting one another like in the attached image. When I run feature to polygon in ArcGIS Pro, I get 7 polygons.
How do I tweak the analysis so I only get 3?
I'm also versed in Python if there is a programmatic way to go about this.


texts: Find nearest boundary of inner feature (Turf.js)


texts:
I have a large shapefile  (~10 Gb) which it output of Atoll (RF Coverage) and I need to perform some post-processing on it.

I began by importing the file in a file geodatabase, the file was imported without problems, but when I run a geoprocessing it raises errors of invalid geometries and the famous 999999. I tried to execute repair geometry but in vain.

I tried to import the shapfile in a PostGIS database using shp2pgsql but during the importing the process stops.

I am stuck on the problem for more than a week, any ideas of how to solve the problem?


texts: print polygon from "the_geom" field in postGIS database

texts:
I have a postGIS database with plots and features (owner, area, address, etc) and I would like to print any selection of that like a card.
I mean I would like to show each record in a page with the bunch of datas and a draw of the polygon that represent the plot . Then I wonder if is possible to read "the geom" field of each record for drawing it.
If not, I would like to know the best method to do that.


texts: How to convert OSM data to a projected Coordinate system

texts:
What is the easiest way to get OSM data into a shapefile projected data format for QGIS?
Sounds simple but I've run into problems with this before.


texts: ArcGIS JavaScript API: Does current extent match initial extent

texts:
In using the JS API, I need to detect if the current extent matches the initial extent. If no match, then I need to display a UI element ("Return to Original Extent").
The initial extent is cached in a data store after the map loads (in my own object mapState.initialExtent). This is properly configured as an instance of an esri extent. I can trigger the map to move to this extent successfully.
I have created a listener on the map (onExtentChange) which calls a function which compares the current extent to the initial extent. It seems the best way to do this is to compare the x/y min/max which is fine and works EXCEPT that BOTH the return from the call to map.extent AND the extent return by the listener are inaccurate.
That is: if I do this:
// SET LISTENER ON MAP LOAD

dojo.connect(map, "onExtentChange", compareExtent)

// AT SOME TIME LATER DO THIS
// this works
map.setExtent(mapState.initialExtent)

function compareExtent(extent, delta, levelChange, lod) {

// FAILS WHEN COMPARING THE FOLLOWING:
 var currentExtent = map.extent;
 var returnedExtent = extent;
 var initialExtent = mapState.initialExtent;

console.log("currentExtent.xmin", currentExtent.xmin, "initialExtent.xmin",        initialExtent.xmin, "extent.xmin", extent.xmin);
console.log("currentExtent.ymin", currentExtent.ymin, "initialExtent.ymin", initialExtent.ymin, "extent.ymin", extent.ymin);

//these do not match
//HERE IS THE CONSOLE READ OUT
currentExtent.xmin -8476373.507203963
initialExtent.xmin -8476287.515547141
extent.xmin -8476373.507203963

currentExtent.ymin 4412584.280379601
initialExtent.ymin 4412589.057693868
extent.ymin 4412584.280379601

}

So, I checked to see if this listener actually returns BEFORE the map is finished but it does not.
Anyone have ideas?

Edit
I am sorry this is not clear. I set the map to the initial extent, change the extent in response to a UI event and then return to the initial extent: the values do not match though the map appear to be correct.
To be clear, on map init I set the extent to initialExtent - using setExtent(initialExtent) where initialExtent is a cached instance of an extent. I then go to a different cached extent and later return to initialExtent. The console shows the returns.
The value currentExtent.xmin -8476373.507203963 is what is returned by the listener.
The value initialExtent is what is initially cached (which is passed as the argument in both  the first and the final setExtent function - - this value is initialExtent.xmin -8476287.515547141.
The following value is what is returned by map.extent after the listener fires after setExtent(initialExtent):extent.xmin -8476373.507203963
I expect these 3 values to be same.


texts: Linear referencing route events for trajectory data

texts:
I am trying to show route events for depths with trajectory data. The depth data has both positive and negative values (starts 50m above sea level and ends km's below the sea level). The below process was followed:

Convert Well Trajectories to Pts using – Feature Vertices to Points (Data Management) Tool
Create Field in Resulting Points Feature Class called Depth. Populate field using geometry calculator and the Z coordinate of the Pt
Use field calculator to convert the negative depths to positive depths.
Create Routes for the well trajectories using – Create Routes (Linear Referencing) Tool.
Calibrate the Routes with the points created earlier using – Calibrate Routes (Linear Referencing) make the measure attribute the positive depth attribute calculated before.
Now reference your event tables against these calibrated routes using - Make Route Event Layer (Linear Referencing) Tool. Note that an event layer is in memory and will need to be exported or sorted to a new feature class. I prefer sort. Note the depths in your reference tables will need to be positive.
Export/ complete a multiple sort to a new Feature Class using – Sort (Data Management) Tool. Sort on Borehole number and depth.
Compute fields for Eastings, Northing and Depth for each point if required and use in whatever application you want.

This process was applied across a number of borehole trajectories in each case the Z-value for borehole event data was distorted for the start of the trajectory, see example below -

Route events layer
Z/M (Event Depth)

-717.9916/139.1
415.6526/160
386.7207/162
186.5/186.5
191/191

Points used to create route
Z values

27.15
-114.45
-176.53
-181.65
-206.77

I assuming the distortion of the events depth is due to the negative lengths (as per link below) used to create route originally.  Is there a possible work around for this that will avoid having to change the height datum?

ESRI Linear Referencing and negative lengths


texts: Expressions to Find data using QGIS Modeler

texts:
Im new to the expressions within the model builder.
I have built a model already using two points as a Station And Backsight, with a line between them.
Also my model computes a length of the line and its azimuth, (both are in their own fields).
I would like to use these fields after the layer has been created to define some parameters in the Wedge Buffer tool.
The tool asks For an azimuth and length but I would like to use the previously calculated layer azimuth and length
Is there an expression to "fetch" that data from within the table?


texts: Runtime error while adding a new field to the existing feature layer using Python

texts:
I have been trying to add a new field to one of my feature layers. And I keep getting an error and I can't figure out what is going wrong.
I am attaching my code below for reference:
When I changed the spatial reference of the 'tracts10_area', it has created a new feature layer called 'tracts10_area_1' in the mentioned workspace.
I have been trying to add a new field called "Populationdesnity" to my existing layer called "tracts10_area_1". And this is the error that is shown:
input_features = ["tracts10_area"]
out_workspace = r"S:\Grad School\FALL 2020\CEE 424 - GIS for Civil Engineers\Homeworks\HW6 PythonEdition\censusGDB\censusGDB\KingCounty_GDB_census.gdb"
out_coordinate_system = arcpy.SpatialReference(4326)
arcpy.BatchProject_management(input_features, out_workspace, out_coordinate_system)

<Result 'S:\\Grad School\\FALL 2020\\CEE 424 - GIS for Civil Engineers\\Homeworks\\HW6 Python Edition\\censusGDB\\censusGDB\\KingCounty_GDB_census.gdb'>

lyrfile = arcpy.mapping.Layer(r"S:\Grad School\FALL 2020\CEE 424 - GIS for Civil Engineers\Homeworks\HW6 Python Edition\censusGDB\censusGDB\KingCounty_GDB_census.gdb\tracts10_area_1")
    arcpy.AddField_management("tracts10_area_1","Populationdensity","LONG",6,"","","pop_den","NULLABLE","NON_REQUIRED")

Runtime error  Traceback (most recent call last):
File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\management.py",
line 3424, in AddField raise e ExecuteError: The value cannot be a joined table
 ERROR 000840: The value is not a Raster Layer.
 ERROR 000840: The value is not a Raster Catalog Layer.
 ERROR 000840: The value is not a Mosaic Layer.


texts: Merging lots of small portions of line to longer adjacent lines using ArcMap?

texts:
I have a hydrologic system split into 100 meter segments which is what I'm looking for, but I also have a lot of very small segments all over my system.

Is it possible to dissolve/merge those small segments automatically within the closest 100m segment next to it for segments smaller than 10m ?


texts: Changing the vertical height from ellipsoid to a local vertical CRS

texts: i have a set of bathymetric data measured in UTM Zone 35 with the ETRS Datum and having the heights from the ellipsoid due to using an RTK base station. Can i transform these raster data sets from the ellipsoid height to Constanta Height (Epsg:5781 http://spatialreference.org/ref/epsg/5781/ ) by using either ArcGis or maybe in gdal? I have tried with ArcGis 10.5 using the project raster tool and choosing the vertical datum there but nothing seems to happen with the heights.


texts: How to calculate the shortest distance for multiple source and targets using PgRouting?

texts:
I have a table with points and a road network which includes driveways (where you enter the highway).

I want to find the point which is closest to a driveway. So I was thinking I can calculate the distance from all the points to all the driveways (maybe within a specific boundary). And from there I can find which point is closest to a driveway. I have joined the points with the network, so I know which ‘source’ I can use from the network. I also know the ‘targets’.

And I can calculate the distance from one point (source-network) to one target using the source and target columns in the network. Using the following query:

    with subset as (
SELECT seq, id1 AS node, id2 AS edge, cost
  FROM pgr_dijkstra(
    'SELECT gid as id, source, target, st_length(geom) as cost FROM roadnetwork,
    473384, 440012, FALSE, False
  ) as di
)
select sum(cost)::decimal (8,2) from subset

But how can I calculate multiple sources to multiple target? So in the query above I want the numbers 473384 and 440012 to be ‘flexible’.

Or is there any other way to calculate this?

Hope someone can help me with this problem


texts:
I have the lastest PostGIS version. When I query:

select ST_HasBBOX(the_geom) from table_geom limit 1

I get this error:
ERROR:  function st_hasbbox(geometry) does not exist

Why this function does not exists


texts: error invalid output extend in zonal statistic

texts: i have a problem with arcmap10.2...i have a dtm(that now has dimension of the pixel 1000) and another raster whit the same resolution with the code of a specific area..what i want to do is calculate for each area the mean altitude..so i know that with the option zonal statistic as a table i can do that..but when i do it i have this error: invalid output extend (error 010092). i tryed the use the default parameters in the environment..but i have the same result..what i should do?


texts: Move node feature and snap to node of other layer

texts: This seems like I'm missing something obvious - I have a layer with node features that I need to update the position of by referencing another node layer.  Should be a simple move and snap-to operation, but there's no snap settings for nodes??? Just vertex and segment?  Is there something that I'm missing?  I was able to do a few with the QAD plugin, but that's been crashing on me and I've got hundreds more to do.


texts: Phenology divergence color map per pixel

texts:
My aim is to create a map displaying the start of spring for vegetation which is a set average value, I want for each pixel to know on which date it is actually trespassing this value. So lets say pixel A is 4 days in advance to the average to trespass this value and is therefore colored in light blue, while pixel 2 is trespassing date of spring two days later than the set average and therefore colored in red.

How can I then create a "divergence colour shp file" which displays a color on each pixel regarding its range or distance to a set value?


texts: ArcGIS API for JavaScript: Draw a Tracker Line from GPS points

texts:
I have gps points which sent from tacker device I need to draw a line in openstreet map ArcGIS with Javascript, How can I do this. I can convert gps points to any format that needs to be (XML, JSON, GPX...) and I have time with lat and lon if needed.

I am new with ArcGIS API for Javascript.

34.673876,68.790897
34.673876,68.790897
34.673876,68.790897
34.673876,68.790897
34.673876,68.790897

I want something like this but ArcGIS its google:

http://jsfiddle.net/pintu31/EPM6A/2/


texts: Counting overlapping polygons in GeoJSON to result in new set of polygons?

texts:
I have 100 polygons in a GeoJSON file that partially overlap each other. I would like to produce 100 new polygons as follows:

The area covered by any of the original 100. This is their union.
The area covered by at least two of the original 100.
...
The area covered by all but one of the original 100.
The area covered by all of the original 100. This is their intersection.

(I'm actually going to only keep 5-10 of these 100 results, but I expect the process for generating all 100 of them is the same as generating a subset of them)

What software / functions can I use to do this?

One approach I've already encountered in another similar project involves rasterizing all of the original polygons, then counting the overlapping pixels and producing new polygons from that data. This approach loses precision earlier than I would like to (I want to maintain near perfect precision up until the end of the process, THEN simplify things for display).

I could do this semi-manually with hundreds of union and intersection operations. This approach is not ideal because I want to re-run this process at intervals in the future as the dataset grows.


texts: Getting from a DEM to points that could be used for a TIN

texts:
I have a DEM, and would like to create a set of points that I can use to make a mesh for a 3d modeling program, such as Maya, Blender, or SketchUp. If I were to make a TIN from a DEM in ArcGIS, ArcGIS would create a greater density of points in areas where the slope changes.

How should I perform a similar process using open source tools? Do you know what algorithm is used for sampling points in a DEM to create a TIN?

I'm looking for a solution that would produce less points than the number of pixels, and in locations where the DEM is more flat, produces less points. The best solution for me would use python or command line tools and result in either a mesh-like geometry or points that could easily be triangulated to create a mesh.

Thanks.

[EDIT: would classifying areas by some combination of slope and aspect be a first step? maybe then I could put points along the borders of these areas?]


texts: Is sen2cor available as a Python library or can it only be used through the command line?

texts: I am looking to correct older Sentinel 2 data to have BOA reflection and would like it to be part of a Python script (I'm working in Windows), but looking online it seems that you can only run sen2cor through the command line. Is there a way to run it directly in a Python script?


texts: Selecting date for torque animation load using Torque?

texts:
I have an animation working using Torque, and the animation is set to pause on load. this automatically loads on the first date of the data series. Is there a way to draw the animation so that it loads on a specified date? My viz also has a datepicker that I would like to be able to set to control the animation so that if I select a date from the datepicker, the animation jumps/redraws and loads to that same date while still having the full timeline of data.

Here's my current code.

Is there a way to set the step in the .on('load') method?

// Draw animation layer using torque
cartodb.createLayer(map, reservoir_storage_layer, options = {https:true, time_slider:true})
  .addTo(map, 1)
  .done(function(layer) {
      layer.on('change:time', function(changes){
          var currentDate = layer.getTime();
          //console.log(currentDate);
          if (changes.step === layer.provider.getSteps() - 1) {
                  layer.pause();
          };
      });
      layer.on('load', function(){
          layer.pause();
      });


texts: Crop OpenStreetMap data retaining long features

texts: When I crop OpenStreetMap data using Osmosis with the --bounding-box option, it seems to only give me nodes and ways that are fully contained within the bounding box. There are several long highways through my region that don't show up at all because no segments lie inside my (relatively small) query area. I don't mind if the long features extend beyond the bounding box, only that all of the features that overlap the bounding box be included. Keep in mind I'm using Osmosis to chop down a large osm file (a cloudmade California file); I'm not running an XAPI query.


texts:
I'm extremely new to the GIS field, and have been tasked with purchasing a handheld GPS to work with a future purchase of ArcGIS for Desktop 10.3.

I need to know if I can upload points (lat,long, and point attributes) from a Trimble Juno 5B to ArcGIS without purchasing Trimble Pathfinder for *.ssf conversion or Trimble GPS Analyst extension for ArcGIS for post-processing?


texts:
I´m trying to find out more on the possibilities with GeoClouding and PostGIS.
I was experimenting with AWS RDS feature with an Free Tier account, and now I´ve been charged a large amount of money. I might misunderstood something here. I´ve been charged for theese RDS-services (which i thought was inlcuded in the 750 free hours per month):

Amazon RDS Service for PostgreSQL ($0.910 per RDS Second Gen Extra Large Multi-AZ Instance hour (or partial hour) running PostgreSQL)

Amazon RDS Service Storage ($0.238 per IOPS-month of Multi-AZ Provisioned IOPS and $0.298 per GB-month of Multi-AZ Provisioned IOPS Storage)

These were tested in 2 regions: EU (Frankfurt) and US West (Oregon) Region.

I´ve only set up an empty  DB-instance,  and now they wanna charge me around 700 dollars.

So the question is: Did I misunderstood/didn´t understand the "contract" well enough or do I have any chance to avoid this? The bill hasn´t been charged from my card yet. It might be my own fault, but it´s jungle of texts and if you´re not from a country with english as your main language, then it´s confusing.


texts: Projection where each pixel in a raster has equal height and width (in km/miles)

texts:
I'm looking to make a raster image where each pixel has an equal width and height in terms of km/miles. My shapefile is in WGS84, so if I convert that straight to a raster, the pixels in the image closer to the equator will be wider (in terms of miles on the surface of the earth) than the pixels closer to the poles. I'd like to avoid this.

Other distortions are OK - I don't care about angles being preserved as much, for example.

Is there a commonly used map projection that supports this?

Edit: Area in question is USA (continental + Alaska, Hawaii, and Puerto Rico), so it's covering a decent portion of the globe.


texts: I have a general question for those who use Esri software: Is there a big difference between the basemaps available for QGIS and those provided free by Esri to its customers? I tried to use the Esri Satellite base map on QGIS and I didn't find a big difference in the ground resolution (houses, streets, etc..), but I guess it depends on the fact that I use another software.


texts:
I'm trying to generate tiles for a region with a surface of 85k km2. The data source are MrSID files (2734 files). I'm using gdal to achieve this task.

It works fine with a small subset of MrSID, but with all the source files it would take more years than I could live :).

My process for generating tiles is the following:

# create the mosaic
gdalwarp -of GTiff -dstalpha -r cubic [MrSID Files] mosaic.tiff`

# generate tiles from mosaic
gdal2tiles.py -r cubic -s EPSG:23030 -z 3-17 -v mosaic.tiff

At first I've chosen a subset of 30 MrSID and it takes 10 hours and it works. Check here.

After that, I've tried with all the MrSID files and I can't reach to create the mosaic, after two days gdalwarp is at 1% processing the first MrSID.

As a workaround, I've tried to avoid the creation of the big mosaic, I've called gdal2tiles for each MrSID, and merging the shared tiles with imagemagic. The issue with this method is that I need a huge quantity of temporal files (30 MrSID files (90MB) give me around 30GB of temporal files).

Any ideas how can achieve this? Is it possible create such a big mosaic? Do you suggest me another strategy?

MrSID files are like this:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: test/807-1-3.sid
Size is 7355, 4796
Coordinate System is `'
Origin = (309892.000000000000000,4291525.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  ICC__Profile=0,0,1,148,65,68,66,69,2,16,0,0,109,110,116,114,71,82,65,89,88,89,90,32,0,8,7,212,0,19,0,2,0,29,0,23,97,99,115,112,77,83,70,84,0,0,0,0,110,111,110,101,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,246,214,0,1,0,0,0,0,211,45,65,68,66,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,99,112,114,116,0,0,0,192,0,0,0,36,100,101,115,99,0,0,0,228,0,0,0,119,119,116,112,116,0,0,1,92,0,0,0,20,98,107,112,116,0,0,1,112,0,0,0,20,107,84,82,67,0,0,1,132,0,0,0,14,116,101,120,116,0,0,0,0,40,99,41,32,50,48,48,52,32,65,100,111,98,101,32,83,121,115,116,101,109,115,32,73,110,99,46,0,100,101,115,99,0,0,0,0,0,0,0,28,69,115,99,97,108,97,32,100,101,32,103,114,105,115,101,115,32,45,32,71,97,109,97,32,50,46,50,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,88,89,90,32,0,0,0,0,0,0,243,81,0,1,0,0,0,1,22,204,88,89,90,32,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,117,114,118,0,0,0,0,0,0,0,1,2,51,0,0
  IMAGE__COMPRESSION_BLOCK_SIZE=512
  IMAGE__COMPRESSION_GAMMA=2.000000
  IMAGE__COMPRESSION_NLEV=8
  IMAGE__COMPRESSION_VERSION=2,0,0
  IMAGE__COMPRESSION_WEIGHT=4.000000
  IMAGE__CREATION_DATE=Sat Jan 12 13:20:09 2002

  IMAGE__INPUT_FILE_SIZE=35277072.000000
  IMAGE__INPUT_FORMAT=TIFF w/ World File
  IMAGE__INPUT_NAME=H:\Imagenes\MrSid\Vuelo_Americano\comp\807_1_3.TIF
  IMAGE__STATISTICS_MAX=255
  IMAGE__STATISTICS_MEAN=231.585752
  IMAGE__STATISTICS_MIN=1
  IMAGE__STATISTICS_STANDARD_DEVIATION=58.039863
  IMAGE__TARGET_COMPRESSION_RATIO=10.000000
  PShop__ImageResources=
  VERSION=MG2
Corner Coordinates:
Upper Left  (  309892.000, 4291525.000)
Lower Left  (  309892.000, 4286729.000)
Upper Right (  317247.000, 4291525.000)
Lower Right (  317247.000, 4286729.000)
Center      (  313569.500, 4289127.000)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Gray
  Minimum=3.000, Maximum=255.000, Mean=232.406, StdDev=53.985
  Overviews: 3678x2398, 1839x1199, 920x600, 460x300, 230x150, 115x75, 58x38, 29x19


texts: Line pattern fill missing when exporting as SVG in QGIS Print Layout

texts:
When exporting from Print Layout as a .svg-file in QGIS the line fill pattern disappears of all polygons. How can I stop this happening?
The problem doesn't occur when exporting as an image file or PDF, however, I would like the map to be in SVG format to edit the layers further in Inkscape.


texts: OpenLayer2 Tween animation callback not working when the browser tab is inactive

texts:
I have data dumping in my db after some interval from resources.The geometry from the data then calling by JS interval function after every ten seconds.using this function to animate my feature http://dev.openlayers.org/releases/OpenLayers-2.13.1/doc/apidocs/files/OpenLayers/Tween-js.html. i which the size of the point change via callback function and represents like animation. All things are working fine when I am on the browser tab. But when i am not on my application tab and during that period if some features draw on my WMS layer then they did not animate normally but stuck like into there begin size. i have tried to console something in mention callback variable. but the behavior i have faced is that all my callbacks kind of stuck or in the waiting mode. When i come back to my browser tab or activate my tab that waiting consoles executes then immediately.
// tween animation
function animatLayer(layer) {
    var animationTween = new OpenLayers.Tween(OpenLayers.Easing.Quad.easeInOut);
    var begin = {
        x: 20
    }
    var end = {
        x: 2
    }
    var callbacks = {
        eachStep: OpenLayers.Function.bind(animate, layer),
        done: function () {
            animationTween = new OpenLayers.Tween(OpenLayers.Easing.Quad.easeInOut);
            begin.x = (begin.x == 20 ? 5 : 20);
            end.x = (end.x == 20 ? 5 : 20);
            console.log('callback')
        }
    }
    animationTween.start(begin, end, 80, {
        callbacks: callbacks
    });
}

function animate(delta) {
    for (var i = this.len; i < this.features.length; i++) {
      this.features[i].style.strokeWidth = delta.x;
      this.features[i].style.pointRadius = delta.x;
    }
    this.redraw();
  }

The call is not happening when the browser tab is inactive so in that case the point got stuck into its max position and does not animate.


texts: Qgis, display asc data file as scatter plot

texts:
I have a digital elevation model in Qgis, and I would like to know how to display the data as points (the X Y grid, with the Z altitude) ? So I can see the density of points in my zone of interest.

The data are an .asc file. And I imported them as a raster layer in Qgis.


texts: Sorting layers by name in one specific group of QGIS layer tree

texts:
I'd like to reorder the layers within ONE specific group of the TOC alphabetically by layer name, but without reordering the toplevel layers or the layers of other groups.
I found the following Python-script:
from collections import OrderedDict
root = QgsProject.instance().layerTreeRoot()
LayerNamesEnumDict=lambda listCh:{listCh[q[0]].name()+str(q[0]):q[1]
                                   for q in enumerate(listCh)}

mLNED = LayerNamesEnumDict(root.children())
mLNEDkeys = OrderedDict(sorted(LayerNamesEnumDict(root.children()).items())).keys()

mLNEDsorted = [mLNED[k].clone() for k in mLNEDkeys]
root.insertChildNodes(0,mLNEDsorted)
for n in mLNED.values():
    root.removeChildNode(n)

The script works principally fine, but it adresses only the top level of the TOC, what is not suitable in my case.
Does anybody has an idea how to modify this script for matching my purposes? Or maybe also a complete alternative?
(I just have some basic knowledge in R, but not in Python, unfortunately)
My conditions in detail:

The top level layers should NOT be sorted, but remain unchanged.
Exclusively the layers of one specific group should be sorted alphabetically, in case of the example only 'group1'  (see example image below)
The layers of other groups ('group2' in the example) should NOT be sorted and also remain  unchanged.

The example is made with a couple of temporary layers. Actually I have to handle a satellite time series with a few hundred multi band raster images, which should be sorted by date.
By the way, I'm aware of the loadthemall-tool. But this would not be useful to me, because I have to add further images with time.

Thanks a lot in advance for any response!
Update 27.5.2021, referring to my comment of 27.5.2021


texts: Linear regression between every 3×3 pixels between two rasters using R

texts: I have two rasters of precipitation and elevation, PPT and ELEV. They have equal number of rows and columns.I want to create a new raster of precipitation lapse rate (GRADIENT) with same number of rows and columns as PPT and ELEV. I want to compute this raster  by performing linear regression between every 3×3 pixels of the two rasters (PPT and ELEV) such that each pixel of the GRADIENT will hold the lapse rate value obtained from linear regression of the corresponding 3×3 pixels in PPT and ELEV that surrounds that pixel.


texts:
I have a dataset of lines with a start and end elevation value in the attribute table. I want to apply this value to the Z coordinate of the start/end vertex respectively to create a 3D line. The lines themselves do not represent a single Z value like a contour; each line represents a sloped downward and thus has a different elevation at its start vertex than at its end vertex.
I did this once in Arc Desktop using 3D Analyst tools by exporting all the start and end  nodes to point shapefiles, converted the 2D points to 3D point using the start/end elevation for each shapefile, assigning all start points with a sequence of "0" and all end points with sequence "1", then merged the 2 shapefiles into one to get a "from to" array of points, finally ran points to lines to get 3D lines.
I need to do this in QGIS because I do not have access to 3D analyst anymore.
Anyone out there done this and willing to share the procedure?


texts: Why can't I find SaveToLayerFile output?

texts:
I have a feature class that has data for the whole month of January with 'date' as one of the attribute categories (labeled ARGH here). I want to grab the information for each day in January and stick it in it's own feature class, so I have 31 separate files. So, as you'll see in the code below, I perform a query for each day, make a feature layer for that day's info and save the layer file to my gdb.

I don't get any error messages when I run it. HOWEVER, I can't find the saved layer files anywhere. I've looked all over. I can't figure out where I've gone wrong.

import arcpy, os
arcpy.env.workspace = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb"
arcpy.env.overwriteOutput = True

Source_File = "AIS_Jan2016_CorrectTIme"

field = "ARGH"
Day_List = [row[0] for row in arcpy.da.SearchCursor(Source_File,"ARGH")]
Day_List1 = [x.encode('UTF8') for x in Day_List]
UniqueDay = set(Day_List1)
UniqueDay_List = list(UniqueDay)
print UniqueDay_List

arcpy.env.workspace = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb"
outputpath = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb"
for day_query in UniqueDay_List:
    out_layer = str(day_query)+ "lyr"
    out_layer_name= "s" + out_layer
    #out_layer_fc= "fc" + out_layer

    #MakeFeatureLayer variables
    in_features = "AIS_Jan2016_CorrectTIme"

    where_clause = ['ARGH = ' + "'" + day_query + "'"]
    #so as an sql-type query this would be ARGH = 'day_query' or something similar
    print where_clause
    workspace = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb"

    try:
        #Execute MakeFeatureLayer
        print in_features
        print out_layer
        print where_clause
        print workspace
        arcpy.MakeFeatureLayer_management(in_features, out_layer, where_clause, workspace)

        #Execute SaveToLayerFile
        arcpy.SaveToLayerFile_management(in_layer, out_layer_name)

    except:
        print arcpy.GetMessages()


texts: Transfering topographic map from ArcMap onto Garmin GPS?

texts:
I have a Garmin 64S GPS unit that I'd like to get a topographic basemap layer onto, so I'm not wandering with a blank screen while getting to my field sites for my research.

I'm using ArcMap 10.3.1, and I have a USGS topographic map that I'd like to download onto my GPS unit.

I've tried using Basecamp and Map Tiler software to get this topo onto my GPS via ArcMap, and I've also tried directly sticking it into a "CustomMaps" folder, but nothing has worked.

My GPS unit, which has the most recent software update, only accepts JPEGs for custom maps, which I've attempted to import to no avail.

I've called Garmin Customer Service, and they wanted me to georeference the .jpg via Google Earth, save it as a .kmz file, and then import it onto the GPS. Again, my GPS doesn't recognize .kmz files. I'm working with a GIS specialist, and we're both scratching our heads currently on how to do this.

Is it possible to take a topo from ArcMap and transfer it onto a GPS unit, or is that impossible?


texts: QGIS 2.2 problem, loading layer into road graph

texts:
I have the road graph plugin installed, but when I try to configure it with the Settings dialogue, no layer names appear in the Layers field. I have tried various shapefiles, but the Settings dialogue doesn't see any of them. Anyone else seen this ? Thanks.

Addition:
I am having the same problem. A few more details regarding my experience of the issue, since not too many are provided above:

-Road Graph plugin version 0.1 with QGIS 2.2.0 32-bit via OSGeo4W install on a 64-bit Windows 7 machine. Tried same software versions on another Win7-64 machine (thinking maybe some setting was tweaked wrong). No difference. I have tried shapefile and spatialite line layers of topologically sound street networks. Tried having just one layer at a time loaded as well as multiple, with no change. Made sure CRS for layer and project were matched and that 'on the fly projection' was disabled. Also de/re activated plugin.

Under Vector>Road Graph>Settings all 3 dropdown menus under the Transportation Layer tab are grayed out/inaccessible with no options shown. All other settings in this dialogue are accessible.

Clicking on the canvas to choose Start and Stop points registers coordinates and shows selected points on the canvas. When 'Calculate' is attempted, "Plugin isn't configured" error (not surprising since layer couldn't be chosen).


texts: Calculating area using dissolve and merge but getting shapefiles with empty attribute tables in QGIS?

texts:
I have a shapefile of a country with the 29 provinces. I also have shapefiles with the areas which are protected from construction (there are 4 categories for protection, within 4 shapefiles. Some of the protected areas belong to more than one category). I have also created a shapefile which is the buffer of the electrical transmission lines.

I would like to work out the area within each province which is maximum 20km from transmission line, minus the protected areas. (This is so I can eventually work out the maximum area available for solar panels.)

I am using QGIS.

My method is like this:

Dissolve the shapefiles for protected areas (as they are broken up by name, etc, but I am not interested in this, only the total area). I have done this by creating a new column ("One") in the attribute table that has the value "1" for each entry and dissolved by this attribute because I cannot figure out a better way. I have then deleted all other attributes since I think they are slowing things up when I am doing other things (maybe this is wrong?)
Merge the 4 dissolved protected areas. With the new layer I dissolve again to get rid of any overlapping features. I am now left with a shapefile of the protected areas which has an empty attributes table.

The next step is where I am getting stuck. This is what I expected to work:

Clip The protected areas layer by the transmission buffer layer. I hoped that this would leave me with a shapefile that has only the protected areas that lie within transmission line buffer area. However this does not work because the resulting layer has an empty attribute layer and therefore cannot be seen.

What am I doing wrong here?

I would eventually like to:
- Create a shapefile of the non-protected area within the transmission buffer area (I expect this would be by using Difference?)
- Split this by province (I have tried to use Join by Attribute but it is not working)


texts: Finding distance along a line (layer 1) where polygon intersects (layer 2) using QGIS

texts:
Layer 1 Line
Layer 2 Polygon

The polygon will intersect the line in 2 places.

E.g.

I'd like to attach the 2 distance values along the line to the polygon layer.


texts: Watershed delineation with meandering river?

texts:
I'm using ArcGIS (Hydrology toolbox) for watershed delineation for a number of pour points, and I'm running into problems with very meandering rivers on extremely flat alluvial land (Amazon basin).

When I run the tool I get extremely small watershed rasters, or larger ones that break off at random points. I have used the exact same DEM, Flow Accumulation and Flow Direction raster files to successfully delineate other watersheds in more mountainous areas.

I'm assuming that the flat relief and meandering river path are to blame, but I'm not sure. The raster resolution is 90m, and I'm using a single point shapefile with the pour point manually placed dead centre on the river. (As I've successfully done previously.) I've tried with a raster pour point but there's no difference.

Does anyone have experience delineating watersheds for this type of river, or any suggestions?