texts: Can't convert geometry from Dgn file to shapefile?

texts:
I am using below codes to try convert geometry from a Dgn file to shapefile. But, I can't exit the Feature reading loop.

Does anyone know why and how to fix it?

Please download the .dgn file from the link: DgnFile

bool Convert()
{
    OGRRegisterAll();
    OGRDataSource       *poDS,  *poSHPDS;
    const char *pszDriverName = "ESRI Shapefile";
    OGRSFDriver *poDriver;
    poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(
                  pszDriverName );

    // Open file DGN
    poDS = OGRSFDriverRegistrar::Open( sourceFile, FALSE );

    if( poDS == NULL )

    {
         cout << "Error opeing DGN file." << endl;
         return false;
    }
    // Create datasource for file shp
  poSHPDS = poDriver->CreateDataSource(destinationPath,NULL);

  if( poSHPDS == NULL )
  {
        cout << "Error creating shape file." << endl;
        return false;
  }
  OGRLayer  *poLayer, *poSHPLayer;

  // Get Layer from DGN file

  poLayer = poDS->GetLayer(0);

  // Create new Layer of SHP file

  poSHPLayer = poSHPDS->CreateLayer (shpFileName, NULL, wkbPolygon, NULL);

  if(poSHPLayer == NULL)

  {

          cout << "Error creating Layer of Shape File." << endl;

          return false;
  }

  OGRFeature *poFeature;

  poLayer->ResetReading();

  **while( (poFeature = poLayer->GetNextFeature()) != NULL ) // Can't exit this loop**

  {

        OGRGeometry *poGeometry;
        poGeometry = poFeature->GetGeometryRef();

        // if feature is a Polygon
        if( poGeometry != NULL && (wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon ) )
        {
                OGRFeature *poSHPFeature;
                //Create new Feature shape file
                poSHPFeature = OGRFeature::CreateFeature (poSHPLayer->GetLayerDefn());
                poSHPFeature->SetGeometry (poGeometry);
                OGRPolygon *poPolygon = (OGRPolygon*) poSHPFeature->GetGeometryRef();

                poLayer->SetSpatialFilter(poGeometry);
                poLayer->ResetReading ();

         }
         OGRFeature::DestroyFeature( poFeature );

    }
    OGRDataSource::DestroyDataSource( poDS );
    OGRDataSource::DestroyDataSource( poSHPDS );
    return true;
}


texts:
While trying to rotate the end caps of a water network in ArcMap, i calculated the rotation of the connected water main using the field calculator, then did a spatial join to transfer the attribute to the cap and used rotation in symbology to implement the rotation. My problem is as shown in the picture below. The end caps that are left of the watermain are correct but the ones on the right are not. It is quite logical if you consider the angle.

The python code used for the angle calculation is:

Expression:
GetAzimuthPolyline( !Shape!)

Code Block:
import math
def GetAzimuthPolyline(shape):
 radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
 degrees = radian * 180 / math.pi
 return degrees

Is there a better way to do this kind of symbology without being affected by the line direction?

Probably something like "perpendicular to line".


texts:
I am using the Landscan grid to calculate the number people within Vietnamese communes. First I masked the data with the GADM boundary. Then I used extract() to calculate the population by commune and saved as a sp file. I converted sp file into a data frame. The code is below.
pop <- raster("w001001.adf")
crop <- crop(pop, gadm)
mask1 <- mask(crop, gadm) # this is so that the next step is a bit quicker
vn <- extract(mask1, gadm, fun=sum, na.rm=T, sp=T)
dta <- as.data.frame(vn)

sum(dta$mask1) # sum of data frame ~ 82 millions
cellStats(mask1, sum) # sum of raster value ~ 95 millions

I think the sum of either the raster or the data frame should give me the entire population of Vietnam and therefore, be identical. However, the two sums are different. The sum of the raster values give me 82 millions (which is most likely correct) and the sum of the data frame yields 95 million people... Does anyone know what I am doing wrong?
Update:
cell <-extract(mask1, gadm, cellnumbers=T, na.rm=T)
> sum(duplicated(cell))
[1] 7
> cell[duplicated(cell)]
[[1]]
    cell value
 1332229 47417

[[2]]
    cell value
 1332228 38578

[[3]]
    cell value
 1332227 41256

[[4]]
    cell value
 1329593 27005

[[5]]
    cell value
 1332233 41888

[[6]]
       cell value
[1,] 248323 29696
[2,] 249202 28011

[[7]]
   cell value
 248324 21301

Update 2:
The Landscan data is available at https://landscan.ornl.gov/landscan-datasets. However, there is an approval process that takes a few days before the data can be downloaded. I'm trying to do the exact same thing with the F162004 DMSP Nightlight value which can downloaded from https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html.
The GADM boundary file can be loaded directly in r
library(raster)
gadm <- getData('GADM', country='VN', level=3) # VN


texts: Formatting SQL expressions for ArcPy with ArcSDE for PostgreSQL feature class?

texts:
I'm trying to use SelectLayerByAttribute to find the last point added to a FC.  I figured I could use the objectid field because it's serial.  I wrote the following:

arcpy.SelectLayerByAttribute_management("controller", "NEW_SELECTION", """ "OBJECTID" = (SELECT MAX("OBJECTID") FROM controller) """)

following the notes on found here

However, I get the following error:

  Runtime error  Traceback (most recent call last):   File "",
  line 1, in    File "c:\program files
  (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6494, in
  SelectLayerByAttribute     raise e ExecuteError: ERROR 000358: Invalid
  expression Failed to execute (SelectLayerByAttribute).

I'm using a Postgres geodatabase and I know all the db's require slightly different formatting requirements.  This one I can't figure out though.


texts: How do I plot data stored in an attribute table (really in a database) in qgis

texts:
I have a csv of UK postcodes from Ordnance Survey that looks like

SE1 0AA,10,531869,179427,E92000001,E19000003,E18000007,E09000028,E05000536

When I import to QGIS using 'Add Delimited Text Layer', I can declare columns 3 and 4 as x and y and everything plots nicely with a new layer created.

I want to store the same data in PostGIS so I can join other data on the postcode and then plot resulting attributes in QGIS. The postcode table, with the same fields, has no geometry and I can't see a way to declare the x and y columns as coordinates, so it is shown in QGIS as just an attribute table.

Can anyone explain what I have missed? I apologise for any mangled terminology - this may be why I haven't found any solution on the web yet.


texts: How to figure out how many ROIs are saved the .roi-file?

texts:
I have some .roi-data with stored regions of interest which I created several years ago.
Now I would like to figure out, how many ROIsfor each class are saved in that .roi file.
Open the .roi in Envi was not successful neither with an editor.
Has anyone an idea how to figure out the number of ROIs in the file??


texts:
IS there a way to split blocks in qgis aside from the split features"

for example i want to select a polygon and have it split as required:

is it better to have the lines in a separate layer or have the polygons themselves split ..

EDIT: what im trying to achieve is not splitting the polygon i just want the outline of being split .. that being with a line or by splitting the polygon itself doesnt make a difference

what i need to is speed the process of getting the splitting outline


texts: Polygonizer-master plugin in for QGIS 2.14

texts: I have tried to download and install the Polygonizer-master plugin but this plugin is "Deprecated". I would like to know if there is another plugin with the same functionality.


texts: Leaflet SVG layer is obscuring the underlying GeoJSON interactivity

texts:
I have SVG on top of GeoJSON in my map. How can a leaflet SVG overlay be placed outside of the leaflet-overlay-pane while keeping GeoJSON in-place within the leaflet-overlay-pane? (For example, how could I put just the SVG into say the leaflet-shadow-pane and then turn pointer events off on the leaflet-shadow-pane  making it possible to interact with the GeoJSON?)

Thus-far, when I set:

.leaflet-overlay-pane svg {
    pointer-events: none;
}

it seems to disable to disable pointer events to both the SVG and GeoJSON, which defeats the purpose.


texts: Having trouble finding out what coordinate system a map is using, and how to convert it

texts:
I'm using Small Area shapefiles in Ireland from the Central Statistics Office (CSO) website. The link to these is here http://census.cso.ie/censusasp/saps/boundaries/ED_SA%20Disclaimer1.htm (Look for the small area ones)

I can't seem to convert the coordinates to latitude longitude, can anyone see what I'm doing wrong?

I've loaded these in using R maptools package using the following code:

cte<- readShapePoly("Shapefile/Census2011_Small_Areas_generalised20m.shp")

But when I check the cte@proj4string it reurns NA. Then I found this document which suggests the projection used is EPSG:3035
http://www.cso.ie/en/media/csoie/census/documents/census2011griddataset/1%20Km%C2%B2%20Grid%20dataset%20User%20Guide%202011.pdf

I then found a site which gives the proj4 string for EPSG:3035. I then tried to set the coorinates as this string.

cte@proj4string <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs" )

Finally I tried to transform the object to latlong coordinates using the following:

cte <- spTransform(cte, CRS("+proj=longlat +datum=WGS84"))

I also tried:

coordinates(cte) <- project(cte, "+proj=longlat +datum=WGS84")

The first transformed it to new coordinates, but they don't seem to be correct. The second one just put everything as infinite.

Can anybody let me know what I'm doing wrong?

Thanks.


texts: Problem inserting feature into GeoServer empty layer via WFS-T

texts:
I'm facing a strange issue with GeoServer (2.4) (on PostgreSQL 9.4/PostGIS 2.1).

The problem is I am not able to insert a feature via WFS-T if the layer is empty; it returns me the following error:

<ows:ExceptionReport
    xmlns:ows="http://www.opengis.net/ows"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://0.0.0.0:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
    <ows:Exception exceptionCode="InvalidParameterValue">
        <ows:ExceptionText>Error performing insert: Error inserting features</ows:ExceptionText>
    </ows:Exception>
</ows:ExceptionReport>

The same WFS-T request works if the layer already has at least one feature.

I also have to note that I'm working with layers based on PostgreSQL updatable views (and I cannot do in another way) using the gt_pk_metadata_table table as reference for PKs.

This is some testing I did:

created new layer based on the view
tried WFS insert transaction (doesn't work)
inserted a feature directly to PostGIS via QGIS (result: the layer has exactly 1 feature)
tried WFS insert transaction (still doesn't work)
in "Edit Layer" page I clicked "Reload feature type"
tried WFS insert transaction (NOW IT WORKS) (result: has exactly 2 features)
removed all features (2) directly via PostGIS (pgAdmin) (result: 0 features)
tried WFS insert transaction (IT STILL WORKS) (result: 1 feature)
in "Edit Layer" page I clicked "Reload feature type"
tried WFS insert transaction (DOESN'T WORK ANY MORE)

During the testing I started a GIT repository in the data_dir to check for file changes but NO file was changed!

I have the feeling it has something to do with new PKs in the updatable view but could be completely wrong.

Can anybody help me figure out what the problem is?

Update:

I changed the logging level to GEOTOOLS_DEVELOPER_LOGGING, as suggested by @user30184, and it looks like the error is due to un-matching srid:

ERROR: new row for relation "table_name" violates check constraint "enforce_srid_the_geom"

I should also note that the WFS-T request is started by OpenLayers using a different projection from the one used in the PostGIS table (EPSG:900913 -> EPSG:3045).

I guess the origin of the problem could be in the fact the layer is based on a PostgreSQL updatable view, which is empty when the layer is empty and I think GeoServer isn't able to figure out the srid of the view thus generating a PostGIS query based on the request srid. Then when the updatable view actually tries to write the record PostGIS realizes the srids aren't matching and fails.

Could this be the reason for the problem? Any suggestion on how to solve this?


texts: does someone know where could I find a description of CRS EPSG:6372. I searched in the EPSG page but seem not to find something that quick..


texts: How to transform any GeoTIFF to Lambert Azimuthal Equal Area with GDAL?

texts:
I am trying to transform a GeoTIFF EPSG:4326 +proj=longlat +datum=WGS84 +no_defs
and a GeoTIFF EPSG:31982 +proj=utm +zone=22 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs both to EPSG:102020 +proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

Command used:

gdal_translate -a_srs "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -of GTiff "source path" "target path"

Results:

The target file for EPSG:4326 seems to be in the right place, but in wrong scale (pixel size unchanged). The target file for EPSG:31982 seems to be in the right scale, but in wrong place.

I already tried the following commands:

gdalwarp -overwrite -s_srs "+proj=longlat +datum=WGS84 +no_defs" -t_srs "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r near -multi -of GTiff "source path" "target path"

and

gdalwarp -overwrite -s_srs "+proj=utm +zone=22 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -t_srs "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r near -multi -of GTiff "source path" "target path"

Both results seems to be in wrong place and wrong angle (rotated).

How to solve this geotransformation with GDAL?


texts:
I am working on Landsat 8 satellite image and want to perform land cover classification for a particular region
Preferred software tool is R
I have formed a RGB image from the satellite data and now want to extract the metadata of that image to work further on it

How can the metadata be extracted?
Can it be done in R or should I use QGIS?


texts: Determining CRS from given point coordinate set

texts:
I have a set of coordinates (X, Y) and I do not know which coordinate system they are actually in.
They are definitely in a projected coordinate system since they are not in decimal degrees. However, I need to determine the CRS and accordingly transform it because they do not overlay other layers.
Below are the sample of coordinates. Perhaps they are in Gauss Krüger CRS. How can I transform it to EPSG:25832?


texts:
I have an ArcSDE 10.0 enabled SQL Server 2008 R2 database on a development server. We are using TeamCity to run continual deployment to this development server, which runs through the following steps:

The production database is restored from a recent backup.
ArcSDE deregistrations are run to remove a layer that isn't used anymore.
Some SQL scripts are executed to make any schema changes that are not in the backup.
Lastly, some ArcSDE registration are run to add new layers.

The deregistration script that runs has no problems:

E:\ESRI\bin\sdelayer.exe -o delete -l jobsite,activityarea -i sde:sqlserver:[server]\[instance] -s [server]\[instance] -D [database] -u [redacted] -p [redacted] -N

ArcSDE 10.0  for SQL Server Build 685 Fri May 14 12:05:43  2010
Layer    Administration Utility
-----------------------------------------------------
Successfully Deleted Layer with Table jobsite.

The registration scripts are both having issues:

Executing E:\ESRI\bin\sdetable.exe  -o create_view -T spatial.vw_seasoncell -t RS_AOI_Grid_Final -c RS_AOI_Grid_Final.shape -i sde:sqlserver:[server]\[instance] -s [server]\[instance] -D [database] -u [user] -p [password]

ArcSDE 10.0  for SQL Server Build 685 Fri May 14 12:05:43  2010
Attribute    Administration Utility
-----------------------------------------------------
    Error: DBMS table not found (-37).
    Error: Unable to create view spatial.vw_seasoncell

And secondly... (note that JobSiteActivityArea is an existing SQL table with a spatial column ActivityArea)

Executing E:\ESRI\bin\sdelayer.exe -o register -l jobsiteactivityarea,activityarea -e p+ -C Id -t GEOGRAPHY -E [extent] -i sde:sqlserver:[server]\[instance] -s [server]\[instance] -D [database] -u [user] -p [password] -G 4283

ArcSDE 10.0  for SQL Server Build 685 Fri May 14 12:05:43  2010
Layer    Administration Utility
-----------------------------------------------------
Error: DBMS table not found (-37).
Error: Unable to create registration for table jobsiteactivityarea.

The real point of interest is if I run these commands remotely (i.e. in my own command prompt - I have ArcSDE on my local machine) targeting that server, everything works correctly. In my mind that would seem to rule out any permission issues with the credentials I'm supplying. This works from a couple of different machines.

I've tried:

Different casing for the tables, as I know they can be temperamental. The JobSiteActivityArea table is pascal cased, but again - it works when it's run remotely.
Adding the schema to the table names. These tables and views are all in the spatial schema, and the credentials that I'm supplying is for the 'spatial' user (i.e. my user is the same name as my schema, because I know it's fussy about that too).

Any suggestions would be greatly appreciated. I'm hoping I'm missing something obvious.


texts: QGis export part of map (box with coordinates) as raster

texts:
I would like to export a map of Switzerland (as part of googlemaps) to use as a background for a graph, but I need to know exactly where each corner of the map corresponds to (wrt spatial coordinates). Is there a plugin to do this in QGIS?

Ideally I would like to be able to manually select a bounding box (while using the OpenLayers underneath) and get a picture exported together with the bounding box coordinates.

Thanks


texts: How to call GDALGetRasterHistogram( ) from R?

texts:
I need to calculate a histogram (pixel-count by each unique value, actually) for more than 256 values , the default number of buckets in gdalinfo.
I understand that GDALGetRasterHistogram can deal with this, but I can't find an option to call/use that function/tool from R o command-line. Any idea?
Here an example:
## Raster with more than 256 values
r <- raster::raster(matrix(1:10000, ncol = 100))
raster::writeRaster(r, 'C:/temp/10000.tif')

## Fast GDAL implementation, but limited to 256 values
currentHist <- capture.output(gdalUtilities::gdalinfo('C:/temp/10000.tif', hist = TRUE ) )
## R option to calulate the values, but very slow in big files
desiredAndVerySlowHist <- table(r[])

Thanks so much


texts: Trying to auto generate the coordinates of a point using QGIS and QField

texts: I am using QGIS 3.10 for windows, and I am trying to create a shapefile such that I can load it into QField and take points in the field, and have those points auto generate their coordinates. I have tried using the $x and $y expressions, but they either are not working or not saving, because every time I try to create a point the Easting and Northing fields are still NULL and when I go back into the Field Calculator the expression is blank. Below are my attribute table and attributes form showing the expression:


texts: with qgis, can you highlight a single feature within currently selected features?

texts:
In Qgis, I have four features that i've selected. They are overlapping polygons and I'd like to cycle through them to see what each feature looks like by clicking on the rows in the attribute table. By default, if you click on a single row, it changes the selection to that feature alone and deselects the other three.

If i remember correctly when this action is performed in ArcGIS desktop, a secondary highlight of that feature is displayed but the original selection remains. If there a way to accomplish this in Qgis?


texts: How to set up a Workspace setting in Geoserver using Curl

texts:
I have successfully created a workspace + store + raster layer in Geoserver using Curl in PHP.

If I go back to Geoserver's UI and enter the Workspace that I have just created, the Settings: 'Enable' and any of the 'Services' are UNchecked (empty). I would like to enable them by using Curl.

If you visit:
http://docs.geoserver.org/latest/en/user/rest/examples/curl.html#

Apparently the only -simplest?- way of creating a workspace in geoserver is by using the XML tags:

<workspace><name>your name</name></workspace>

If I try to add new tags to the XML expression above, such as for instance

<enabled>True</enabled>

it gives you an error.

[500] No such field org.geoserver.catalog.impl.WorkspaceInfoImpl.enabled

So, Have anyone successfully configured those options of the Workspace by using Curl? For instance enabling WMS Services for that workspace?

By the way this post does not help too much:

Selecting service (WMS / WFS) through GeoServer REST API?


texts: Why does the File Geodatabase that I created not appear to any other users on the network?

texts:
I created a File Geodatabase to store files that I use in maps for my coworkers.  The folder is in a location where all users should have access and permissions.

Strangely enough, after I created the gdb, I'm the only one that can see or access it.  It appears in my Catalog view, but not in windows explorer. The map document that accesses the fcs in the gdb works fine for me but shows broken links for my colleagues. Any reason why this might be happening??


texts: Matching data attributes between records of different layers

texts:
I'd like to know how I can match data attributes between records of different layers, in which there are the same fields, in one shot.

For example, i'd like to transfer data from "test_1" to "test_2" quickly, especially when there are lots of fields. Any suggestion?


texts: Is there a way to read the columns in a shapefile / layer without converting it first to an MDB file?

texts:
My shapefile has a number of fields, like roadname, and street number.  Is there a way to get at these without converting the shapefile first into a featureclass (mdb)?  Sometimes the conversion is a time consuming operation.

I'm using ArcEngine 10 C# with VS2010


texts: How to merge multiple shapefile via ArcEngine 9.3 Java?

texts:
I try to merge the some shapefiles into one by using ArcEngine java platform.
  But it seem like that parameters of my merge class is wrong ,i google it a long time and     got nothing. Pls figure my problem out .The Error is :Invalid value type for parameter inputs.

public static  IGeoProcessorResult mergeTool(String shpPath,String[] shpNames,  String outFile)
{
    //Initialize the geoprocessorResult object
    //IGeoProcessorResult vGeoPR = null;
    // Initialize the geoprocessor.
    GeoProcessor vGeoP = null;
    IGeoProcessorResult vGeoPR = null;
    try {
        vGeoP = new GeoProcessor();
        IWorkspaceFactory vWSpaceFactory =new ShapefileWorkspaceFactory();
        //IFeatureClass vFeaCls = vFeatureWSpace.openFeatureClass(ShapeName);
       IWorkspace vWorkSpace;
       IFeatureWorkspace vFeatureWSpace = null;
       vWorkSpace = vWSpaceFactory.openFromFile(shpPath, 0);
    vFeatureWSpace = (IFeatureWorkspace)vWorkSpace;
        IFeatureClass vFeaCls = null;

        ArrayList<IFeatureClass> vInArray = new ArrayList<IFeatureClass>();
        for(int i=0;i<shpNames.length;i++)
        {

            vFeaCls= vFeatureWSpace.openFeatureClass(shpNames[i]);
            vInArray.add(vFeaCls);
        }
        //vGeoP.addMessage("I will overwrite your file if name exits!");
      // vGeoP.setEnvironmentValue("workSpace", vFeatureWSpace);
        vGeoP.setOverwriteOutput(true);
         Merge  vMerge = new Merge(vInArray,outFile);
            vGeoPR=vGeoP.execute(vMerge,null);
            return vGeoPR;


texts: How do I calculate Manhattan distance with PostGIS?

texts:
I am using the ST_Distance function to calculate the distance between two geometries (a train station and a building) Since I know that all buildings and all train stations are in Chicago, which has an excellent/complete street grid, I would like to use Manhattan (or taxicab) distance.

The generic formula for this is the difference in X plus the difference in Y, so Abs(X1-X2)+Abs(Y1-Y2).

What PostgreSQL query would make this work?


texts:
How can I display a UTM Grid as a layer with Leaflet?
Search a way to embed a OpenStreetMap map with Leaflet in a website.
Particularity. I want to display a UTM Grid as a Layer.
In addition, lower right UTM coordinate where the mouse is.

Example with UTM grid and UTM Number:
http://geo.dianacht.de/bayluftbilder/utm.php

Solution?
https://github.com/kartena/Proj4Leaflet

Does anyone an example?


texts: How to Clip a DEM and retain DEM format?

texts:
I have a DEM file that I downloaded from the internet.  I am trying to clip it to reduce the file size.  So far I have tried two different methods, but neither of them produces a clipped file in DEM format (which I need to pass on to an AutoCAD user in DEM format).

I have tried the Clip Tool under Data Management -> Raster -> Raster Processing.  I used another image file for my Output Extent.  The problem is that if I use DEM as the extension under "Output Raster Dataset" it will not run the tool.  I am able to run this tool if I choose TIF, but I need it as a DEM.
I have also tried right-clicking the DEM and choosing Data -> Export Data.  I drew a graphic rectangle that I wanted to clip to, and used Selected Graphics (Clipping) for my extent.  With this method, DEM is not an option in the Format drop down menu.  I can use this method to export to a TIF, but again I am looking for a DEM.

My question is similar to this one, but the answers here don't address how to keep the DEM format:
Clipping raster using vector polygon in ArcGIS Desktop without having Spatial Analyst?

Is DEM a proprietary format (like SID) that I can't export without a special license or program?  Or am I missing something else?


texts: Is there a way to make multiple adjacencies explicit using the Polygon Neighbors tool?

texts:
I am analyzing a shapefile of precincts (polygons) in ArcMap 10.6.1 and used the Polygon Neighbors Tool to generate a table denoting the neighbor relationship between polygons. I am using the field PSN (Precinct Serial Number) to uniquely identify polygons. The neighbor list appears in the following format, with LENGTH being the length of the shared perimeter in meters and NODE_COUNT indicating whether the neighbor touches only at a point. The neighbor list goes both ways; so if 533 and 534 are neighbors, there are instances of both precincts showing up in the src_PSN field as seen below.

In my map, I have several (I'm unsure of how many) multi-adjacent polygons; in other words, polygons that touch their neighbors in two places. For example, here is one instance in my map. The two highlighted polygons, 533 and 538, have this multi-adjacency.

Is there a way to make multi-adjacencies show up as two separate records?

In the example above, polygon 533 would have 11 neighbor records rather than 10. I need the multi-adjacency to be explicit in my table in order to run a bit of code I've been handed to run an analysis on the map.


texts: How to use Collect Values and Tabulate Areas in ModelBuilder?

texts:
I have problems with the Tabulate Area tool.

The polygons in my shapefile overlap and because of this ArcMap does not calculate proper values.

I try to solve this using ModelBuilder where it is possible to do the Tabulate Area for each polygon separately but at the end I would like to combine rows back to one table.

I thought that the solution may be the model only tool called Collect Values?

Can someone help me how to add rows to one table or to save each Tabulate Area result under different names?

I hope that the question is understandable. If not, I will try to rewrite it.

Here is also a picture.


texts: Fixed inverted Buffer in QGIS for 100 meters is not working

texts:
In QGIS I am trying to run a buffer on a polygon layer for 100 meters, it is based on a tutorial so it is supposed to work but it doesn't work.  Someone in another forum suggested that I merge the polygons first but it actually didn't make a buffer at all after that.

I am using WGS84 UTM Zone 33S, the tutorial link is here http://docs.qgis.org/2.14/en/docs/training_manual/answers/answers.html#complete-analysis-2


texts: How to delete by attribute from shapefile in SAGA command

texts:
I have lot of shapefiles to erase polygons with same attribute form them. It is usually an easy task when do manually, but how can I use SAGA-GIS command?
Any help?


texts: Match housenumbers on buildings (special case of point-in-polygon-test)

texts:
Task with example
I'm working with geodata (country-size) from openstreetmap. Buildings are often polygons without housenumbers and a single point with the housenumber is placed within the polygon of the building. Buildings may have multiple housenumbers.
I want to match the housenumbers to the polygons of the buildings.

Simple solution
Foreach housenumber perform a point-in-polygon-test with each building-polygon.
Problem
Way too slow for about 50,000,000 buildings and 10,000,000 address-points.
Idea
Build and index for the building-polygons to accelerate the search for the surrounding polygon for each housenumber-point.
Question
What index or strategy would you recommend for this polygon-structure? The polygons never overlap and the area is sparsly covered.

This is a duplicate from stackoverflow.com. It was recommended to post the question here.


texts: How to display a Shapefile on a browser? -offline-

texts:
Plot points on a Shapefile like this: https://github.com/gka/pyshpgeocode
and then convert it to HTML and display it on a browser?

Is this practical method to display the map on a browser?
Any cons of it?

Is there any other way of displaying a Shapefile on a browser?


texts: How to clip raster and split by individual polygons in QGIS

texts:
I have two layers. The first one is a vector layer with 6 polygons. The other one is a raster layer. I need to extract specific surface information for my polygons (average slope value). The information is available in my raster layer « Alma-ShipshawPente ».

I want to clip my raster image so it fits the extent of my vector layer and is split by individual polygon features. How should I proceed?

I looked at a few tutorials to help me through. All of them included programming. I am not familiar with programming and, therefore, I have not succeeded yet. Are there other ways to proceed ?

I tried the option Raster > Extraction > Clipper. Here's the result.

It is a beginning, but the raster is not split by polygon features. It is clipped according to the extent of the entire layer.


texts:
What fields can the CalculateField_management expression parameter contain? I mean the parameters can be the fields what are in the output table, or can they be a field of another table too?

There are two tables:

lot of fields and a field, what needed to be placed to other table.
The joined table. It contains the field what comes from the previous table and all of the records has unique value in that certain field.

Now the task is to go through on the table 2's certain field and compute other field in that table 2, but I need to search for values in table 1. If a record's "certain field's" value is the same as table 2's "certain field" value, then I need the other field datas from table 1 to compute table 2's other fields.

Does the expression parameters need to be from only one table or can I refer to other table fields too?

I hope it's understandable.


texts:
I have a map in the form of .bil & .hdr using WGS84. I'd like to convert this into another map using a 2D projection while keeping the lat/long information. What I'd really like with this new map is to have a consistent distance/pixel. With the WGS84, the distance/pixel changed from top/bottom and left/right. The map size is about 100 miles.  Here are my question.

Is this possible to do this to a level of accuracy of +/- 100 ft.
What would be a good 2D projection? I have seen people recommending not using mercator but I dont know what my other options are.
What is the command to transform this using gdal. I think this is it but the # of options were a bit overwhelming and I wasnt sure if I was doing this right.

gdalwarp -t_srs '+proj=merc +datum=WGS84' input_wgs84.bil out_merc.bil


texts: Is it possible to manually change Network position of Network Analyst Stop without moving the stop?

texts:
Network Analyst does not always choose the street I need.
When I get a new .mxd file, the Route territories have been defined and the Network Analyst Stops placed centered at the nearest network street (within the county parcel shape) or where the position has been marked from a previous visit.
When I Solve, the solution does not always take me to the correct road.
The Network Analyst Stop graphically represents the target, but sometimes the house and driveway are not on the same street. Sometimes the Property and Parcel are landlocked with no network street associated.
Without moving the stop from the house (target) to the driveway, how can I manually change the position on the network?
This is difficult in rural areas where the roads and driveways are thousands of feet apart, but also in downtown areas. A house at 101 5th Avenue will be routed around the corner to Baker Street.  Since there is a sidewalk easement on the front of the parcel, the Stop is closer to the side street than to its address.


texts:
I've been attempting to download a usable basemap for USGS Topography for quite some time. Going through the USGS procedure to acquire the map has not worked to date.

Quickmapservices and Openlayer plugin have only given me a black layer that does not work.

I'm also trying to download a shapefile for FEMA Floodplain Zones. Its been a very frustrating process and even after calling and going through the steps, nothing seems to work.


texts: Deploy Google App Engine app for Google Earth Engine

texts:
Basically, an app from Google App Engine makes a request to Earth Engine (server to server basis) then renders this image to a web page.

In my case, it requested SRTM data:

    mapid = ee.Image('USGS/SRTMGL1_003').clip(bound).getMapId({'min': 0, 'max': 1000})

I was able to replicate demo described in Google EE demo apps in my localhost and online as well: https://service00-201702.appspot.com/

However, as you see the DEM is in black and white and without legend.
Do you know where to config its color ramp and legend which shows elevation range with color ?

This is github link for this app: github.com/google/earthengine-api


texts:
Have a JavaScript application trying to push data to a FeatureServer on ArcGIS Server 10.1.

I get this error occasionally:
{Error code '400': 'Unable to complete operation.'
Lock request conflicts with an established lock [dbo.DEFAULT]}


texts: Seperate building-layer into buildings within villages and remote buildings using QGIS

texts: im trying to seperate a building-layer into two layers. One layer should show the buildings within villages and the other one the remote buildings. My first approach was to buffer the buildings (and dissolve) so i can seperate it by the area of the bufferareas. Big areas are villages, small ones are the remote buildings. It failed because the dissolve function did dissolve the attributes too. The second approach was to create a grid and change the polygons of the buildings with a pointlayer so i can count the points inside a cell of the grid. Unfortunately I could not find a function that works (im beginner with GIS-Softwares). Do any of you have an idea how I can solve this task?


texts: GeoServer greyscale image with alpha not working

texts:
With GeoServer 2.16.4 I have the following problem:
A greyscale GeoTIFF renders fine in layer_preview which is defined (gdalinfo)
Size is 11994, 5020
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 31N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 31N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 0°E to 6°E - by country"],
        BBOX[0,0,84,6]],
    ID["EPSG",32631]]
Data axis to CRS axis mapping: 1,2
Origin = (520026.485833066864870,4034738.386493506375700)
Pixel Size = (20.000255836264788,-19.999740259740243)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.9.dev with GDAL 3.0.4
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  520026.486, 4034738.386) (  3d13'24.60"E, 36d27'27.69"N)
Lower Left  (  520026.486, 3934339.690) (  3d13'15.45"E, 35d33' 9.04"N)
Upper Right (  759909.554, 4034738.386) (  5d53'56.32"E, 36d25'21.68"N)
Lower Right (  759909.554, 3934339.690) (  5d51'57.91"E, 35d31' 7.12"N)
Center      (  639968.020, 3984539.038) (  4d33'10.25"E, 35d59'42.98"N)
Band 1 Block=11994x16 Type=Byte, ColorInterp=Gray
  NoData Value=0

The same image with a set alpha Channel does not render
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  520026.486, 4034738.386) (  3d13'24.60"E, 36d27'27.69"N)
Lower Left  (  520026.486, 3934339.690) (  3d13'15.45"E, 35d33' 9.04"N)
Upper Right (  759909.554, 4034738.386) (  5d53'56.32"E, 36d25'21.68"N)
Lower Right (  759909.554, 3934339.690) (  5d51'57.91"E, 35d31' 7.12"N)
Center      (  639968.020, 3984539.038) (  4d33'10.25"E, 35d59'42.98"N)
Band 1 Block=11994x16 Type=Byte, ColorInterp=Gray
  NoData Value=0
Band 2 Block=11994x16 Type=Byte, ColorInterp=Alpha
  NoData Value=0

My guess was that it's the colorspace. However converting to RGAB or YCBCR did not solve it either. I would be thankful if one could advise me how to render the image with an added alpha channel in GeoServer.
Find here the full stacktrace: https://hastebin.com/comusetune.log


texts: Has the "Warp (Reproject)" tool changed in QGIS 2.12?

texts:
I am writing a QGis plugin, for QGis 2.10, which imports processing 2.10.3, that needs:

testproc = processing.runalg('gdalogr:warpreproject',
                          thisName, # input
                          thisCRS, # source crs
                          thisNewCRS, # dest srs
                          '0', # no data, <parameterString>
                          0, # target resolution: 0=unchanged
                          0, # method: 0, as we are only clipping
                          0, # output raster type
                          2, # compression
                          None, # jpeg compression
                          None, # zlevel
                          None, # predictor
                          None, # tiled
                          None, # bigtiff
                          None, # TFW
                          extraParam, # extra
                          thisOut)

extraParam is a string I used to clip my raster with a shapefile, using -clip option.
It works for me with processing 2.10.3. I need to give this plugin to someone using processing 2.12.2. with QGis 2.12.0 (since this version of QGis requires a more recent version of the plugin). It fails: I think that the input for this command changed.

how I can detect the version of processing (after doing import processing). I would like to adapt the behaviour to the version, as the long term release uses 2.10.3
what changed in the new version, in term of inputs?


texts: Loss of vector line accuracy when exporting from QGIS

texts:
I'm new to QGIS and having a little trouble understanding why my vector lines in QGIS look sharp and accurate but when I export and then import into MapBox, the line accuracy is reduced.
Vector end points shift a little and don't meet correctly anymore. There seems to be some vector simplification happening during export?
I've tried multiple export formats and they all provide the same result. Any ideas?


texts:
I'm also having problems with running gdal_merge.py within another python script and I cannot run it. I already checked a similar problem here. My example is as follows:

First I set my workspace and locations:

sys.path.append('C:/Python27/ArcGIS10.2/Scripts/')
import gdal_merge as gm workspace="D:/Satellitendaten/rapideye/img/testregion/cannyedge/out/"
os.chdir(workspace)

Then I run:

sys.argv = ['o', 'out.tif', 'allre1.tif', 'allre10.tif', 'allre11.tif', 'allre12.tif', 'allre13.tif', 'allre14.tif', 'allre15.tif', 'allre16.tif', 'allre17.tif', 'allre18.tif', 'allre2.tif', 'allre3.tif', 'allre4.tif', 'allre5.tif', 'allre6.tif', 'allre7.tif', 'allre8.tif', 'allre9.tif']

I get:

  0...10...20...30...40...50...60...70...80...90...100 - done.

But also ERROR 4:

  `out.tif' does not exist in the file system, and is not recognised as a supported dataset name.

So the output is not created. Maybe someone can help me out with the issue to merge different tif files to a single tif?


texts: How can you convert between geographic data from and to 3D data in PostGIS?

texts:
Is there was a way to convert a geography value to a 3D geometry value (in ECEF, for example), and vice versa? I know how to convert between geography using a 2D projection, but not a 3D reference frame.

I could probably create a function to do this myself (using this and this, for example), but is there already a function to do so, (hopefully more general or efficient than I could write)?


texts:
I have some issues with Quantum Gis 1.8. When I try to use shapes to grid (sextante-saga tool) I get the following error (comes up when I try to define the new grid file name):

Fehler bei der Ausführung von Python-Code:

Traceback (most recent call last):
  File "C:/Users/Max/.qgis//python/plugins\sextante\gui\OutputSelectionPanel.py", line 94, in saveToFile
    settings.setValue("/SextanteQGIS/LastOutputPath", os.path.dirname(str(filename)))
UnicodeEncodeError: 'ascii' codec can't encode character u'\xfc' in position 19: ordinal not in range(128)

Python-Version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS-Version:
1.8.0-Lisboa Lisboa, f440b60

Python-Pfad: ['C:/Users/Max/.qgis//python/plugins\sextante_taudem', 'C:/Users/Max/.qgis//python/plugins\sextantelwgeomprovider', 'C:/Users/Max/.qgis//python/plugins\sextanteexampleprovider', 'C:/Users/Max/.qgis//python/plugins\sextante', 'C:/Users/Max/.qgis//python/plugins\permaclim', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/Max/.qgis//python', 'C:/Users/Max/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\PROGRA~2\Quantum GIS Lisboa\bin\python27.zip', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\DLLs', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\plat-win', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\bin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\python\plugins\fTools\tools', 'C:/Tempordner_GIS/J\xc3\xbchnde']


texts: Does the experimental QGIS2leaf plugin not work on mac OS?

texts:
I am trying to get the experimental QGIS2leaf (Leaflet web maps) plugin to work on my macbook air (2012, Mavericks OS), but when I open the outputted index.html file, no interactive image shows up in my web browser (I have tried Chrome, Safari, and Firefox).

Is this a common issue with the QGIS2leaf plugin? Is there a work around?


texts: Create shapefile based on geojson data with geotools (Java)

texts:
I'm using Geotools 14.0 for Java and I'm trying to generate shapefile based on data from GeoJSON. Here is my code:

try {
        String url = "geojson.json";
        File geojson = new File(url);

        File shpFile = new File("test.shp");
        ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();

        Map<String, Serializable> params = new HashMap<String, Serializable>();
        params.put("url", shpFile.toURI().toURL());
        params.put("create spatial index", Boolean.TRUE);

        ShapefileDataStore shpDataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);

        InputStream in = new FileInputStream(geojson);
        int decimals = 15;
        GeometryJSON gjson = new GeometryJSON(decimals);
        FeatureJSON fjson = new FeatureJSON(gjson);

        FeatureCollection<SimpleFeatureType, SimpleFeature> fc = fjson.readFeatureCollection(in);

        SimpleFeatureType type = (SimpleFeatureType) fc.getSchema();

        shpDataStore.createSchema(type);
        Transaction transaction = new DefaultTransaction("create");

        String typeName = shpDataStore.getTypeNames()[0];

        SimpleFeatureSource featureSource = shpDataStore.getFeatureSource(typeName);

        if (featureSource instanceof FeatureStore) {
            SimpleFeatureStore featureStore = (SimpleFeatureStore) featureSource;

            featureStore.setTransaction(transaction);
            try {

                featureStore.addFeatures(fc);
                for(Object c : featureStore.getFeatures().toArray()){
                    System.out.println("C: " + ((SimpleFeatureImpl)c).toString());
                }
                transaction.commit();
            } catch (Exception ex) {
                ex.printStackTrace();
                transaction.rollback();

            } finally {
                transaction.close();
            }
        } else {
            System.out.println(typeName + " does not support read/write access");
        }
    } catch (Exception e) {
        return "Exception " + e.getLocalizedMessage();
    }

And after this, i get test.shp file but it does not contain any geometry data. If I convert it back to geojson I get this:

{   "type": "Feature",   "properties": {   },   "id": "test.1" }

And this is my original geojson file (geojson.json):

{
  "type": "FeatureCollection",
  "features": [
    {
      "properties": {
        "version": 0
      },
      "id": "1",
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              1868732.467515988,
              6594375.3042187244
            ],
            [
              1868732.467515988,
              5165920.1196253505
            ],
            [
              3297187.6521093622,
              5165920.1196253505
            ],
            [
              3297187.6521093622,
              6594375.3042187244
            ],
            [
              1868732.467515988,
              6594375.3042187244
            ]
          ]
        ]
      }
    }
  ]       }

Is there any way to do this in a proper way?


texts: Python error in polygon to raster conversion?

texts:
I want to convert polygon to raster. I converted the dbf file ( "cntymam") to table view to make use in the conversion. I am getting a error that the parameters are not called and also with the grid name. I am not able figure out the error.

Errors

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000161: The length of the grid name must not exceed 13 characters
ERROR 000800: The value is not a member of CELL_CENTER | MAXIMUM_AREA | MAXIMUM_COMBINED_AREA.
Failed to execute (PolygonToRaster).

Script

# Make layer and table view

arcpy.MakeFeatureLayer_management(FLcounties, "flcnties")
arcpy.MakeTableView_management(cntymam2, "cntymam")

j=0

for species in mamList:
    if j < 10 :
        species_str = str(species)

        # polygon to raster

        mamcnty_out = output_workSpace + "cnty_" + species_str

        Value_field = "cntymam2." + species_str

        print species_str, mamcnty_out, Value_field

        arcpy.PolygonToRaster_conversion("flcnties",Value_field, mamcnty_out, "CELL_CENTRE", "NONE", "90")

        print j
        j= j + 1


texts: Grouping by with multipolygon geometry with sf (using R)

texts:
I have a custom shapefile I put together. It works fine when I plot it all at once. But I'd like to group by certain variables to plot specific shaped regions. For example:

+-------------+--------+-------+
|   county    | region | sales |
+-------------+--------+-------+
| washoe      |      1 |     5 |
| carson city |      1 |    10 |
| clark       |      2 |    15 |
| harmon      |      2 |    20 |
+-------------+--------+-------+

If I run:

leaflet() %>% addTiles %>% addPolygons(data=df)

That will plot four separate polygons based on where those four are located.

But let's say I wanted to create polygons based on regions. So theoretically, the output would contain just two Polygons. One that is the combined shape of montgomery and stevens, and another combining moses and harmon county, and I also want to sum up the total sales.

x <- df %>% group_by(region) %>% summarise(totals=sum(ones))

It will give me this error:

Error in CPL_geos_union(st_geometry(x), by_feature) :
  Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point -119.97497995771261 39.521600907266631 at -119.97497995771261 39.521600907266631

In terms of information on the geometry variable:

df$geometry

Geometry set for 2106 features
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -115.8968 ymin: 35.00184 xmax: -114.0428 ymax: 36.85366
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Do I need to convert it to something else?


texts:
First sorry if I am incorrect in the words, I am using a translator to expose my doubt.
I have a requirement where I have a worldwide road base and my goal is to calculate the road network by country.
I know that it is impossible to do this calculation using only half a dozen "datums and projections", but I would like to know if there is a list of which EPSG to use by country or even by regions.
Here in Brazil we have the Brazilian Institute of Geography and Statistics - IBGE which suggests which EPSG to use for calculating the extension of linear elements (EPSG: 5880).
The big question is:
Is there any documentation where you can find a list of which EPSG to use for each region?
If we do not have it ready I believe that we are already late in creating a large repository where the best EPSG to use is listed.


texts: weighted sum overlay of raster layers (Wind resource map and GHI Map)

texts:
I want to generate an output raster map which indicates the areas of wind and solar energy potentials. My study area is India.
I have created GHI (Global Horizontal Irradiance) raster map (by interpolating point GHI values available). Similarly, I created wind speed map,by interpolating MERRA satellite data.
Now for making Wind-Solar hybrid map, if I add the two raster layers by "Raster Calculator" of Spatial Analyst Tool in ArcGIS 10.1, then this is the right way of doing this?


texts: I have a question which has been asked before. I've been using QField very successfully on an android tablet, but have bought a new one and can't open a project, although it is in the same place as in the earlier tablet - Android/data/ch.opengis.qfield/files. On the tablet I cannot find the Android directory. The new tablet is a Samsung S2 8 inch with Android 7.0


texts: How to remap a raster to its min / max values in PostGIS?

texts:
I have a table of 16-bit signed (!) integer heightmaps (SRTM tiles) loaded into PostGIS. This website showed me how to export a raster as a PNG to a file:

COPY (SELECT encode(ST_AsPNG(raster), ‘hex’) AS png FROM table_with_raster) TO ‘/tmp/myimage.hex’;

xxd -p -r /tmp/sf_colors.hex > /tmp/sf_colors.png

However, because my image data is in signed integers, ST_AsPng(rast) fails because it's signed integers, but ST_AsPng() needs values from 0 to 256 (8-bit, unsigned). So I want to reclassify the image, so that 0 = minimum value and 256 = maximum value.

For that reason, I created a new table and calculated the statistics of the image with ST_SummaryStatsAgg:

CREATE TABLE heightmap_data_stats as  (
    SELECT rast_id, (stats).min, (stats).max, rast
    FROM (
        SELECT ST_SummaryStatsAgg(rast, 1, true) AS stats, rast_id, rast
        FROM heightmap_data GROUP BY rast_id
    ) AS foo
);

Now I have something like this:

id      min     max   rast
-----------------------------
0       85      1232  [data]
1       37      677   [data]
3       95      625   [data]
2       173     977   [data]

Now I tried:

SELECT encode(ST_AsPng(
    ST_Reclass(rast, 1, 'min-max:0-256', '8BUI', 0)
), 'hex') FROM extra_data.heightmap_data_stats;

... in order to map the image data from 85..1232 -> 0..256 for example. However, it doesn't work - I don't know the syntax to tell PostGIS to input the actual "min" and "max" values in the ST_Reclass argument, so I only get an error "invalid reclassarg" and the original raster is returned.

How do I write a query, that, for each raster in the table, reclassifies it to an 0-256 range?


texts: Projection of raster to sinusoidal projection

texts:
I have my definitive raster shown below. This raster has the projection as mentioned here

class       : RasterLayer
dimensions  : 7276, 8330, 60609080  (nrow, ncol, ncell)
resolution  : 463.3127, 463.3127  (x, y)
extent      : 6508617, 10368012, 751956.5, 4123020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : <address path of where this raster is stored>
names       : <definitive raster>
values      : 0, 578  (min, max)

I have another land use land cover raster (each pixel value being categorical) that is in WGS1984 coordinate system i.e. no projection. I want project this raster to the same projection as the definitive raster to do more analysis.

I did so in R as follows

lulc_reproj<-projectRaster(lulc_, res=res(defintive raster),crs=crs(definitive raster), method='ngb')

It did complete with no errors and I could continue with my analysis. However, I viewed the projected lulc raster and the definitive raster in Arc and see that two rasters do not align. For example, below is a snapshot of the water class in the reprojected lulc raster (light blue) which is way off from the water in the definitive raster (which is NoData and displayed in red)

The same is the case when viewed in QGIS. Like I mentioned earlier, no errors come up in R and the projection strings match additional to the resolution.

I also tried projecting the lulc raster in ArcMap as I have read that Arc uses different transformation algorithms than gdal, on which rgdal and raster packages are based. However, I get the exact same results i.e. no alignment

So why are the rasters not aligned even if they have the same projections and resolutions? How do I fix it?


texts:
I'm wondering if there is a workaround to this.

I want to symbolize a feature class categorically by unique values.  The field that I'm doing this on has about 5 unique values.  I also have a coded value domain on this field with about 3000 values.  When I try and symbolize the field, it creates a symbol for every value in the domain, not just for the values in the feature class.

Is there something I'm missing that will limit the symbology to only the unique values that are in the feature class?


texts: How to enable Latin names for places in some places instead of local ones in Carto Services?

texts:
Some places on the map is labeled with Cyrillic names, but I need only English/Latin names of places on the map, however sometimes there are only local names.
In general I need transliteration to English?
How can I implement this?

P.S.: I have spotted this issue on Belorussian and partly on Russian places.


texts:
I'm trying to create a map with CartoDB to show the spanish population between...at a municipality level.

I can do everything perfectly, but I have just one problem: I would like to know if is there any possibility to customize the "search" feature, connecting it to the geolocalized names that I have in my table. In the way that if I put "Guadalajara" -that I have in my table- it goes directly to Guadalajara, Spain (the one in my table) and not to Guadalajara, México. And, in the same way, If I put "Madrid" right now it goes to the city, just a bit more north than where my point is.

Hope I've been clear. Maybe my solution is to add a filter feature to a map, but I still don't know how to do it.

Thanks!


texts: Do I have to convert a line shapefile to points to store it in another database?

texts: This is the surprisingly the first time I've come across this.  If I wanted to export a line shapefile from ArcGIS to another database (say Access) would I first need to convert the vertices to points and add the x y columns?  To then reuse the data in a GIS I would need to convert back to lines.  I feel sure I must be missing something.


texts: Leaflet layers on different zoom levels. How?

texts:
I've added a polygon layer and a point layer to my leaflet map.

I would like the points to appear after a certain zoom level or when a "click" event occurs e.g. (polygons.on("click", function....)

I have looked at different questions here and also this example

But I am still unable to understand what to do in my code.

Here is my code:


texts: Iterating over a layer and make MBR for every feature with QGIS

texts: I have a layer with about 20000 polygons in it. The make extent layer command just draws the MBR (Minimum Bounding Box) for the whole layer. I want to make MBR for each of them separately. Is there any way to iterate over the layer and make MBR for each of them?


texts: Listing features and their respective map grid in QGIS

texts: I am trying to list map features and the Reference Grid they fall into. The grid is setup in the Print Composer as a 50mmx50mm grid over a 400mmx400mm map with Cells Labeled as 1A, 1B, etc.. I didn't make the Grid as a spatial layer. Is there a way to list the features and the cell they can be found in referencing the Grid directly or do I have to recreate it as a spatial layer? If it has to be a spatial layer then is there an easy way to export the composer grid to a spatial format?


texts: How to convert an ASCII file (with lat,lon coordinates) to a shapefile using OGR?

texts:
I have an ASCII file with the following columns : ID, longitude, latitude, value

Is it possible to use OGR (ogr2ogr or ogrinfo) in command line in order to open the ASCII file and convert it to a shapefile (in WGS84) using the point coordinates (X field=longtitude, Y field=latitude) ?

What I want to do is similar with what QGIS does when opening an ascii file with "Delimited Text" and defining the geometry (point coordinates).

I want to do that for a high number of files so I want to know if this is possible in OGR.

I finally tried the following:

ogr2ogr -f "ESRI Shapefile" test.shp input.csv -oo X_POSSIBLE_NAMES=lon -oo Y_POSSIBLE_NAMES=lat -oo KEEP_GEOM_COLUMNS=NO

but I get an error
Unable to open datasource test.shp with the following drivers


texts: Custom PROJ CRS for shifted, rotated, and skewed axes?

texts:
It is possible to create a custom PROJ definition for a cartesian CRS with shifted, rotated, and skewed axes? For example:

Where The x-axis is formed by the line AC, the y-axis by the line AB, and the origin at A. I'm thinking something like a typical UTM system, but the origin is shifted from the typical origin, rotated so axes don't align with cardinal directions, and skewed so the axes are not at right angles.


texts: Making WFS in Leaflet app hoverable and clickable?

texts:
I have made a WFS on Geoserver for a simple grid shapefile to server to a web mapping app. Within the web mapping app, I want the grid to be hoverable and clickable. Specifically, I want to be able to highlight a particular grid cell when the mouse hovers over that cell and keep the rest of the grid hidden. I have my wfs parameters defined as such and am receiving the json/geojson correctyl in my web app:

 var url = 'http://localhost:8080/geoserver/ows';
 var defaultParameters = {
   service: 'WFS',
   version: '1.1.0',
   request: 'GetFeature',
   typeName: 'FGCMS:1_Deg_Grid_polygon',
   //maxFeatures: 200,
   outputFormat: 'application/json',
   format_options: 'callback:getJson',
   SrsName : 'EPSG:4326'
 };

I want to add this functionality (hover, click, etc...) in the front end in my web app code but I have not found a specific Leaflet plugin that might work (I am not seeing anything specifically in the API) and I do not know how to just add to the json data configuration.

In the code below, I have my ajax request with the grid json configured to fillOpacity: 0. Could I just add something like hoverable: true?

var WFSLayer = null;
var ajax = $.ajax({
  url: URL,
  dataType: 'json',
  jsonpCallback: 'getJson',
  //async: false,
  error: function() { console.log('Oh snap!'); },
  success: function(response) {
     WFSLayer = L.geoJson(response, {
        style: function (feature) {
            return {
                //fillColor: 'FFFFFF',
                fillOpacity: 0
            };
        },


texts: Retrieve attributes from a GeoJSON feature in openlayers 3

texts:
I'm working on a webmap using OpenLayers 3 and vector layers, geoJSON. I've managed to render the layers and add click interactions:

        select = new ol.interaction.Select({

            condition: ol.events.condition.mouseclicked,
            style: myStyle,
        });

        var map = new ol.Map({
            interactions:
            ol.interaction.defaults().extend([select]),

The interaction works fine, it sets a new style to the selected feature (borders of the selected country etc).

What I want is to retrieve the selected feature's id from the geoJSON file and use that to know what country is selected. Depending on what country is selected a network (lines) will be drawn from that country to other countries, retrieved from a DB.

Is there a way of doing this? I've tried the following (without proper formatting):

select = new ol.interaction.Select({

            condition: ol.events.condition.mouseclicked,
            style: myStyle,
        });
var features = select.getFeatures();
var feature = features.item(0);
var id = feature.getId();

This gives me an id that is undefined even though the feature has an id in the geoJSON file.

Is there another way of doing this? I've searched around the web and this forum without finding an answer.

This is a snippet of my geoJSON file, all features have the same properties:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "sovereignt": "Afghanistan", "id": 1 },
"geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 61.210817091725744


texts: Mapping filename into fields of geodatabase using ArcGIS ModelBuilder?

texts:
I'm using ModelBuilder to append a number of shapefiles to a single geodatabase. The shapefiles are located in a folder SHP and the model is designed to iterate through SHP, appending each .shp to the geodatabase. The model enters the name of the file into one of the fields of the geodatabase. Each .shp has very limited information attached to it (field names: FID, Shape, ID).

However, the geodatabase has many fields (Year, Name, Project, Type, etc.), some of which I would like to populate using parts of the filename; the filenames have a specific naming structure, for example: PRJ_YEAR_TYPE.shp, where PRJ is replaced with the name of the project related to the file, and similarly for YEAR and TYPE.

So far, I've been able to add the name of the shapefile to one of the fields of the geodatabase by:

Adding a Field after I bring in my feature class (.shp) and
Calculating the Field (named using one of the preexisting fields in the geodatabase), before appending the file to the geodatabase.

However, I'm not able to add mutliple fields. I've tried to duplicate the processes I used to place the filename into the geodatabase name field, but it doesn't seem to work.

How can I add the inforamtion going into the other fields (i.e. (Year, Name, Project, Type, etc.) making reference to the filename?

I am new to ModelBuilder.


texts: How to remove duplicate labels using ArcPy in ArcGIS Pro

texts: I have a script that uses AddDataFromPath() to add a layer from ArcGIS Online that is regularly updated to a map in ArcGIS Pro. My script then updates and displays the labels using listLabelClasses(). Everything runs fine except for there are duplicate labels. My question is, how do I remove the duplicate labels using ArcPy? I don't see it in the ArcPy documentation so I'm wondering if there is a work around for this. Is it possible to somehow import those settings from a template layer?


texts:
I'm writing a plugin for QGIS (3.10) in Python. I would like to follow the execution of the script line by line to understand some bugs. I installed Remote Debug plugin and tried to connect it with :

PyDev (Eclipse). I configured all as explained here: How to set up Eclipse to develop QGIS scripts or plugins. But executing the plugin in QGIS still has no effect in Eclipse.
Winpdb-reborn, but a window is opening in winpdpb showing the code of Remote Debug plugin and highlighting this line:

def statusBar(self):
    return self.iface.mainWindow().statusBar()

Eric6, and the connection seems to be OK in Eric6, but QGIS freezes as I click on "Connect" in Remote Debug.
First Aid plugin seems not to be compatible with QGIS 3.

I don't know in what direction to go now. Does someone have a working solution to debug plugin in QGIS 3?


texts: How to generalize SpatialPolygons and conserving topology in R?

texts:
I would like to generalize a SpatialPolygonsDataFrame  in R to make plotting more efficient. The polygons I have result in huge pdf plots that are slow to load and hard to integrate into a document.

I tried to use the gSimplify function from the rgeos package, but unfortunately it does not conserve the topology of the polygons.

Is there an alternative that I can use?


texts: QGIS how to handle circular but non-closed lines

texts:
I use QGIS 2.18.3 and I have a shapefile with around 30,000 lines and a lot of invalid geometries. It is imported from a CAD project. I managed to solve most of the invalidities with the break function of the v.clean tool from the GRASSGIS toolset.

The remaining errors, around 180, all have the same nature: they used to be closed lines in the CAD project (the outlines of stones) but are now open. So the two endpoints have exactly the same coordinates, therefore they overlap and the geometry is invalid.

Solutions I've tried:

Used the snap function of the v.clean tool - the tool runs without problems but there are still errors afterwards. I've put the threshold on 0.005 m. In advanced parameters I've put the v.in.ogr snap tolerance to 0.005 m as well and the v.in.ogr min area to 0.0
Used "lines to polygons" from the vector menu and converted the polygons back to lines again. All polygons were valid but when I convert them back to lines using "polygons to lines" the errors are still there.

Did I do something wrong or do you have any idea how to solve it in an automated way?


texts:
Ok since my previous problem of overlaying WFS onto WMS is still unsuccessful, I though was able to overlay the WMS instead with this code,

{layers: 'iPlant:AllQuebecSpecies', transparent: true}

So now I wanted to know can I use the openlayers get-feature example to obtain info from the WMS or does it work only for WFS?


texts: Dissolving/connecting discontinuous Sewerage network using ArcGIS Desktop?

texts:
I am using ArcMap/ArcGIS Desktop 10.3.

I have a sewerage network that is cut at specific locations and I would like to have a continuous network while retaining the different segment because they have different id's.

Someone have any ideas?

0
down vote
accept

The Extend Line Work, but some time it doesn't work correctly like the next example:

Befor Extend line:

After the Extend Line :

Why it didn't take the line just beside-it and it went looking for a line further?


texts: how to run batch process for vector shapefiles

texts:
I am kind of new at this.

I ran a gis operation using python, where by I performed euclidean distance for specific state road, river and towns data and reclassified them showing good and bad areas and it worked.

However when I try to run batch processing for two state fils in the same folder, it runs but no output. any ideas how i can batch process vector files before reclassifying the raster output

import arcpy
from arcpy import env
from arcpy.sa import *
import os

# Check out the ArcGIS Spatial Analyst extension license
if arcpy.CheckExtension("Spatial") == "Available":
   arcpy.CheckOutExtension("Spatial")
else:
 # raise a custom exception
    raise LicenseError

if arcpy.CheckExtension("3D") == "Available":
   arcpy.CheckOutExtension("3D")
else:
 # raise a custom exception
   raise LicenseError
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")
else:
 # raise a custom exception
   raise LicenseError

if arcpy.CheckExtension("3D") == "Available":
    arcpy.CheckOutExtension("3D")
else:
 # raise a custom exception
    raise LicenseError

env.workspace = "D:\\states\\ritchie"
root_folder = "D:\\states\\ritchie"
env.overwriteOutput=1

river = "C:\\states\\ritchie\\river\\"
try:
    fcs = arcpy.ListFeatureClasses()
    for fc in fcs:
        outEucDistance = EucDistance(os.path.join(river,fc), 2000, 50)
        outEucDistance.save("D:\\states\\ritchie\\river\\")
except:
print arcpy.GetMessages()


texts:
I'm working in ArcGIS 10.4.1 and trying to write a script (in Python) for making a polygon shapefile from latitude/longitude.

txtFile=arcpy.GetParameterAsText(1)
shpOUTpolygon=arcpy.GetParameterAsText(2)
shpOUTpoint=arcpy.GetParameterAsText(3)

en=os.path.dirname(shpOUTpolygon)

env.workspace = en
FCpolygon=os.path.basename(shpOUTpolygon)

FCpoint=os.path.basename(shpOUTpoint)

mxd=arcpy.mapping.MapDocument('current')
df=arcpy.mapping.ListDataFrames(mxd)
dfr=df[0]
sp=dfr.spatialReference

arcpy.CreateFeatureclass_management(en, FCpolygon, "polygon","#","#","#",sp)
cursor =arcpy.da.InsertCursor (FCpolygon, ["SHAPE@"])
array=arcpy.Array()
point=arcpy.Point()
for line in fileinput.input(txtFile):
    point.ID, point.X, point.Y=line.split()
    array.add(point)
polygon = arcpy.Polygon(array)
cursor.insertRow([polygon])
fileinput.close()
del cursor

For beginning, I test script on XY coordinates. It works fine.
First, I choose some projected coordinate system manually (in Data Frame Properties of Map), then I run my script and created shapefile  is correct.

But when I input longitude and latitude - no errors while script is working, but something wrong with polygon's geometry in output shapefile.

Where should I set a right spatialReference of polygon?


texts: Display a list field in a custom form with QGIS and Qt Designer

texts:
I'm building a form with Qt Designer and I want to show attributes of my table when I interrogate an object. I design this form with a .ui file.

With integers and strings fields I'm using QLineEdit and don't have any issues. But I can't display attributes with a list field, QGIS crash.

I've tried QListView, QTableView, QListWidget and QTableWidget, same result.

How can I display a list field? Which class of QObject I have to use?

I'm using QGIS 3.10 (LTR) and Qt Designer 5.11.2

EDIT

Here are my form settings. A the level of "creation" field (list) I want to see my multiples values, like in a drop-down menu where the user can scroll.

Here are my QGIS settings. As you can see, "creation" field is a list.

How can I display this kind of list? (form generate automatically)


texts: Searching closest region from Google Maps coordinates

texts:
Here is my problem. I'm mapping a region in Google Maps after that I'm saving the coordinates in PostGIS database as geographic polygon type.

ST_GeographyFromText('SRID=4326;POLYGON((..))')

The goal of my script is to find the closest region to a given point selected from Google Maps.

Here is my solution which is not working very well.

So what I'm doing is to break the region polygon and to find the closest point to my desire point. But the accuracy is not very precise. Some times when my point is closer to a given region it selects a region which is pretty far from it.

SELECT
        ST_Distance (
            ST_GeographyFromText('SRID=4326;POINT(" . $point->longitude . " " . $point->latitude . ")'),
            geom::geography
        ) AS distance
    FROM
        ST_DumpPoints (gis_coords::geometry) --gis_coords is the region column.
    ORDER BY
        distance ASC
    LIMIT 1
) AS distance

When I try to render my region column with QGIS and OpenLayers Plugin the result is pretty disturbing.

QGIS Settings

I'm setting the region layer CSR to WGS 84 (EPSG:4326) and for OpenLayers WGS 84/Pseudo Mercantor (EPSG: 3857)


texts:
I have a large ESRI file geodatabase, that I loaded into QGIS with no problems. If I look at the attributes tab in the properties dialog, the fields are as they should be.

One of the fields is supposed to have small integers. If I try to symbol by this attribute, it gives me four levels: 10005, 48255, and 152528, but almost all of the features are missing the attribute. I cannot open the attribute table without QGIS crashing. Here's a screenshot of the file with the supposed symbology.

If I save the .gdb as a shapefile, then the attributes are properly coded and assigned. Here's the shapefile symbology.

It seems as though there is a problem with how QGIS is handling the attributes table, but that filedb is all right because the shapefile export is fine. I would file a bug report, but I'm wondering if the community has a known workaround or things I could check.

session info

OS X 10.10 Yosemite
QGIS 2.8.0
GDAL/OGR1.11.2


texts: I can't see my WMS layer on Cesium Globe scene

texts:
After insert a WMS layer in my Cesium, I can see in my browser that the images are comming from Geoserver but they are not visible in the globe.

Here is the code ( commented lines are because I'm testing ):

var west = -45;
var south = -24;
var east = -42;
var north = -22;

var rectangle = Cesium.Rectangle.fromDegrees(west, south, east, north);

Cesium.Camera.DEFAULT_VIEW_FACTOR = 0;
Cesium.Camera.DEFAULT_VIEW_RECTANGLE = rectangle;

var imageryProviderCHM = new Cesium.WebMapServiceImageryProvider({
    url : 'http://10.5.115.180:36213/geoserver/wms',
    layers: 'CHM:ElementosProntos',
    parameters : {
        transparent : true,
        tiled : true,
        format : 'image/png'
    }
});

var viewer = new Cesium.Viewer('cesiumContainer',{
    //sceneMode : Cesium.SceneMode.COLUMBUS_VIEW,
    //imageryProvider : false,
    //timeline : false,
});

//viewer.imageryLayers.removeAll();
viewer.imageryLayers.addImageryProvider( imageryProviderCHM  );

var terrainProvider = new Cesium.CesiumTerrainProvider({
    url : '/tilesets/test',
    credit : 'Dr. Who'
});

//var widget = new Cesium.CesiumWidget('cesiumContainer');
//widget.imageryLayers.addImageryProvider( imageryProviderCHM );
//widget.imageryLayers.addImageryProvider( provider );

//viewer.scene.terrainProvider = terrainProvider;

And here we can see the tiles comming from my geoserver and the area they must be:

What I'm doing wrong? Maybe layer order?

EDIT
A clue: My browser is throwing this in console: An error occurred in "N": Failed to obtain image tile X: 777 Y: 640 Level: XX.


texts: Geoserver installed successfully but WFS doesn't work

texts:
I installed GeoServer 2.1 RC1. It installed successfully, but when i click on WFS in services category, i got the following message.

Oops, something went wrong...
Sorry, something unexpected happened on the server.

Here's an error report you can include in a JIRA bug report about this issue:

org.apache.wicket.WicketRuntimeException: Can't instantiate page using constructor public org.geoserver.wfs.web.WFSAdminPage()
        at org.apache.wicket.session.DefaultPageFactory.createPage(DefaultPageFactory.java:212)
        at org.apache.wicket.session.DefaultPageFactory.newPage(DefaultPageFactory.java:57)
        at org.apache.wicket.request.target.component.BookmarkablePageRequestTarget.newPage(BookmarkablePageRequestTarget.java:298)
        at org.apache.wicket.request.target.component.BookmarkablePageRequestTarget.getPage(BookmarkablePageRequestTarget.java:320)
        at org.apache.wicket.request.target.component.BookmarkablePageRequestTarget.processEvents(BookmarkablePageRequestTarget.java:234)
        at org.apache.wicket.request.AbstractRequestCycleProcessor.processEvents(AbstractRequestCycleProcessor.java:92)
        at org.apache.wicket.RequestCycle.processEventsAndRespond(RequestCycle.java:1250)

What is wrong ?
How can i fix this problem ?


texts: GeoTools WFS Query sortBy not working for GeoServer WFS

texts:
I am using GeoTools to query a GeoServer WFS using the following code:

    Map connectionParameters = new HashMap();
    connectionParameters.put("WFSDataStoreFactory:GET_CAPABILITIES_URL", String.format("%s&request=GetCapabilities", wfsUrl));
    connectionParameters.put("WFSDataStoreFactory:TIMEOUT", connectionTimeout*1000);
    DataStore data = DataStoreFinder.getDataStore(connectionParameters);
    SimpleFeatureType schema = data.getSchema("airports");
    FeatureSource<SimpleFeatureType, SimpleFeature> source = data.getFeatureSource( typeName );
    String geomName = schema.getGeometryDescriptor().getLocalName();
    Envelope bbox = new Envelope( -100.0, -70, 25, 40 );

    FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2(GeoTools.getDefaultHints());
    Object polygon = JTS.toGeometry(bbox);
    Intersects filter = ff.intersects( ff.property( geomName ), ff.literal( polygon ) );

    Query query = new DefaultQuery( typeName, filter, new String[]{ geomName } );
    query.setSortBy(new SortBy[]{ff.sort("country", SortOrder.ASCENDING)});
    FeatureCollection<SimpleFeatureType, SimpleFeature> features = source.getFeatures( query );

The query returns features, but sorting never works.  A toString() on the query returns something like:

  feature type: airports    filter: [[ geom intersects POLYGON
  ((-117.8395870053188 31.58066038426718, -115.5957867039628
  31.58066038426718, -115.5957867039628 33.27328332002949, -117.8395870053188 33.27328332002949, -117.8395870053188 31.58066038426718)) ]]    [properties:  ALL ]    [sort by: country ASCENDING]

which indicates that sorting is part of the query.

The following direct WFS query works:

http://thing:8080/geoserver/ix/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=airports&srsName=epsg:4326&bbox=-117.8395870053188,31.58066038426718,-115.5957867039628,33.27328332002949&sortBy=country

Any suggestions?


texts: I was looking for a comparitive statement of the accuracies achived by various GNSS systems such as QNSS, COMPASS, GLONASS, GPS..etc. however i find only discriptions of only GPS everywhere


texts: Precision loss when transforming coordinates systems in GDAL (vector layers)

texts:
Anybody knows the order of magnitude of the error introduced when transforming coordinates using GDAL?

More specifically, I would like to know if significant error is introduced when the transformation involves a change in Datum like, say, from NAD83 to WGS84.


texts:
I would like to make a table (it is enough to display it in the layout) in ArcGIS that would have the attached field (column) structure.
Is it possible and how? I can use ArcPy, preferably.

Edit: Something like MultiIndex for columns
https://jakevdp.github.io/PythonDataScienceHandbook/03.05-hierarchical-indexing.html


texts:
We have a positioning system where we can measure our range without 0.1m from multiple beacons.  The placement of the beacons is about on a 1-2km box.     So we could get anywhere from 2-4 ranges from the beacons.  The beacons are surveyed in place to an accuracy of about 1m.

We found the PHP code for the range case: Trilateration using 3 latitude and longitude points, and 3 distances  and we're testing this in our "C" program w/ double precision math.

Questions:  We're looking to get the most accurate position projection in LAT,LONG over a short area say max range up to 2km.   What's the best aproach to get the least errors in the solution math?  How big will the position errors be in meters?

Is there a code example of finding the two solutions in lat,long is we only have 2 ranges?


texts: Detect lines not connected to the main network R

texts:
How could I autodetect lines not touching the principal net in R? In this example I would like to remove the top right line.

I already know how to select and delete lines by their lengths but that will not work as I need those little lines that work as connectors between points and net.
This is what I'm tying so far: try to detect the intersections between the layer of lines with himself and get rid of those that just have one intersect, what I guess means that they are just intersecting with themselves.

l <- readOGR("./Layers/Lines.gpkg")
st_intersects(st_as_sf(l),st_as_sf(l))

It would be nice any solution.


texts: Using output parameters as input parameters in Python Script Tools?

texts:
I created a script tool in ArcMap 10.5. How do I set an output parameter from a tool in the script as an input parameter in another tool at a later point in the same script?  Like:

buffout = arcpy.GetParameterAsText(5) # Feature class output of buffer analysis-  THIS IS THE VARIABLE I GET AS OUTPUT THAT I WILL LATER USE AS INPUT
buff_dist = arcpy.GetParameterAsText(6) # Buffer distance parameter

# Run the Buffer Analysis tool on stream feature
arcpy.Buffer_analysis(outfc, buffout, buff_dist)

try:
    # Add Parameters
    outrast = arcpy.GetParameterAsText(7) # Output raster dataset clipped to stream buffer polygons
    # Run the clip tool on hillshade using "buffout" variable as output extent
    arcpy.Clip_management(outpath, "", outrast, buffout, "", "ClippingGeometry", 'NO_MAINTAIN_EXTENT')
   # Add a success message
   arcpy.AddMessage("Clipping of hillshade to buffer output successful!")
except:
   # Report any error messages that the Buffer tool might have generated
   arcpy.AddMessage(arcpy.GetMessage(0))

I want "buffout", the output from arcpy.Buffer_analysis() step, to be used as an input parameter in arcpy.Clip_analysis().

Do I define buffout as arcpy.SetParameterAsText()?

If so, how and where do I use arcpy.SetParameterAsText() in my script to accomplish my purpose?

I don't understand that function very well.

Also, I'm assuming I must set buffout as a "Derived" parameter type in the script tool properties after using arcpy.SetParameterAsText() correctly?


texts: Installing a Minimal qgis, without extensions (GRASS, SAGA)?

texts:
Is it somehow possible to install QGIS without the GRASS and SAGA extensions?

I just want to have a small software package to create and export some simple maps and my hdd is quite full already.

I really won't need these extensions but the software package of qgis.com includes all.


texts: How to export data from sql database in the decimal format in QGIS?

texts:
Currently I am working with QGIS and fetching data from SQL server and export it in the format of shapefile in QGIS. While importing the data from sql server attribute table showing data with geography having decimal value.

But after exporting in the format of shapefile data in the attribute table not showing the decimal value.

according to article as follows

decimal degree data erase itself between edit session in qgis

in the description they mentioned to add new column and make decimal places precision as 5 or more than that, but attribute table showing for data imported from sql doesn't have  option to make decimal places and add the precision as 5 or more than that. If we check this option in the attribute table of shapefile we can found that. But it not open properly. It got crashed while editing the column. Another question in front of me is how to copy all huge data to new column?

I wanted to know:

How data got automatically converted in to integer type after exporting it from database
Does any new way to add decimal precision for adding new column in sql database imported attribute table and copy all data in the form of decimal places?
Does any tools, plug in to convert data type from float to float or any mapping type option for data type available?


texts: Memory issues projecting a large raster mosaic of USGS DEMs in R

texts:
I am trying to project a large raster to UTM (meters), but I keep crashing my laptop (16GB) due to a memory shortage. It seems like using projectRaster tends to deplete the machine's memory, and some solutions have been suggested in questions like this one. However, even using the suggested solution my laptop keeps crashing. The raster in question is a mosaic that I create from 4 DEM tiles downloaded from USGS with a resolution of 1/3 arc second (~10m) that correspond to the city of Milwaukee, WI. This mosaic is 396.8 Mb. The process goes as follows and can be run in any machine to download the rasters:
library(raster)

rasterOptions(memfrac = 0.3) # the solution suggested in several questions to limit access to memory by the raster package.

import the four rasters via their individual link.

r1 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w089/USGS_1_n43w089.tif")
r2 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w088/USGS_1_n44w088.tif")
r3 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w088/USGS_1_n43w088.tif")
r4 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w089/USGS_1_n44w089.tif")

# create a list with the four rasters and merge them into one single raster layer
rl <- list(r1, r2, r3 , r4)
mosaic <- do.call(merge, rl)

rm(r1,r2,r3,r4, rl) # remove unnecessary files in memory to save space

# project the mosaic layer to UTM zone 16N (m) -- THIS IS WHERE THE COMPUTER CRASHES
mosaic <- projectRaster(mosaic, crs = "+proj=utm +zone=16 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

In this other question, rasterOptions is used to set chunksize = 1e+04 and maxmemory to 1e+0.6. I tried doing this and my machine didnt complete the line mosaic <- do.call(merge, rl) in the entire night.
I have considered the option of projecting each raster (r1, r2, r3, r4) separately before creaing the mosaic, but I am weary of the issues that this may cause at the time of merging them together, and I would very much rather find a memory efficient way to handle the reprojection than just avoiding it with a by-pass.


texts: How can I use different config.yaml files for application using mapfish print in geoserver?

texts:
System: Ubuntu 10.04 server / Tomcat6 / geoserver2.1.2 with integrated GeoWebCache

I have various webmapping applications running which utilize the mapfish print module.Some are in tomcat6/webapps some are in my apache2 DocRoot. I need to have various mapfish print layouts to show different logos, legends, scalebars etc depending on which webmapping application is using it....but I only have one config.yaml file to style the pdf layout.

Qustion: How can Have I various versions of config.yaml and link them to the different webapps?

I´d be grateful for any replies,

yours,

Rob


texts: Forking and adding to FeatureLayer using ArcGIS Online

texts:
I am trying to find ways to join existing polygon feature layer with new data.  For example, Esri has a living atlas reference layer for US state boundaries: https://www.arcgis.com/home/item.html?id=99fd67933e754a1181cc755146be21ca
This layer also has some population data for each state.  I want to create a new polygon feature layer that is essentially a fork of this one.  I then want to add some data - in this case, voter turnout data for each state.  That data comes in excel format here (source http://www.electproject.org). Once I have a copy of the states polygon feature layer with the new election data, I can create some nice visualizations.  I am trying to create a 2020 version of this feature layer: US Voter Turnout 2016.  But I am having a hard time figuring out how to create and serve a layer from my own account that is a fork of that initial layer, and then add to it.
I have tried going into the Content section of arcgis online, hitting 'create', and creating a feature layer from an existing URL, using the arcgis server URL of the states boundaries item.  When I go through the form and hit 'Done', I get this error:
Unable to add feature service definition. The operation failed because an index or statistics with name 'insights_fl_tauCNbQYU1nWzn8S' already exists on table 'user_72472.lklklkj_USA_STATES_GENERALIZED'.
I have also tried going into the original boundaries layer, and exporting the data to an excel file.  Then I can add my new data as new columns.  Then when trying to create a new layer, I go to create a new Polygon FeatureLayer.  I hit 'Update Data', and then 'Append Data to Layer'.  I try to upload the excel file, but it gives the error: No compatible layer was found in the uploaded file. Please check your file and try again.
I feel like this should be relatively simple to do but I'm struggling to find a way, or any documentation on how to do this.


texts: Calculating the distance between a point and a virtual line of two lat/lngs

texts:
Please refer to the example and corresponding image.

I would like to achieve the following: provide two locations (lat/lng), which are shown below as A and B. From this, a virtual line would be drawn and then the distance between this line and C would be calculated (in any measurement).

I have achieved this currently in Google Maps API v3 but would want to also be able to perform this behind the scenes in my language of choice. Any tips/ideas would be greatly appreciated!


texts: GWR with dependent and explanatory variables in different files

texts:
I have a point shapefile (let's say A) with pollution measures (let's say Y). My hypothesis is that pollution at theses locations is primarily driven by nearby industrial units and not by activities at the locations in A. I have the locations of the factories in another shapefile (let's say B) along with a host of explanatory variables (let's say X1, X2, X3).

I want to model the pollution measure Y as a function of X1, X2, X3, and the distance of the industrial unites from the points of measurements. I was planning on using a Geographically Weighted Regression (GWR), but am not certain whether it is the right choice.

In ArcMap for example, the GWR tool expects one shapefile as input with both the dependent and explanatory variables.Thus, the assumption is that Y is driven by events at that point along with some influence of nearby locations. I looked into GWR in R also and the assumption is similar since it expects only one argument with data. I know that in R I can force feed the explanatory variables and dependent variable to be from different files, like this:

gwr.sel(A$Y ~ B$X1 + B$X2 + B$X3 + B$dist, ...)

But I am not sure if this makes sense, primarily for 2 reasons:

My dependant and explanatory variables are not in the same file. So, I may be violating a basic assumption of GWR.
The dist variable only has the distance from the nearest industrial unit to the measurement location. Thus, I am loosing the effect of the industrial plants which are close to the measured location, but not the closest.

Thus, my question is what other alternative do I have to GWR that will help me model this phenomenon?


texts: QGIS: Delete information from rasters and shapes that are out of country borders

texts:
I am trying to make topographic maps for some countries and show only those countries. I already created reliefs from SRTM-data and imported lake and river data from Natural Earth. Now I want to delete all raster and shape information that is out of the country (in my example down there, Peru). The country is another shapefile from Natural Earth.

Long story short: How do I delete all the stuff out of Peru?


texts: How to make empty parts of aerial imagery transparent?

texts:
I have this problem in QGIS which is very similar to this one: Making empty parts of aerial imagery transparent

I want to connect all these different rasters, so that they form one big aerial image. I already found how to remove the big white borders, so everything looks much better. However, there still remain smaller dots at the borders of the different rasters. Anyone knows how to remove these?

First image:

Second image:


texts: Automating NetCDF data update on GeoServer

texts:
I have published a NetCDF file (received on 2018-01-02, hourly data over the following 7 days) to GeoServer as WMS with Time enabled.

I was able to preview the layer with OpenLayers with the "time" attribute attached to the URL to the GetMap command.

http://{GEOSERVER_DOMAIN}:8080/geoserver/fdi/wms?service=WMS&time=2018-01-02T01:00:00.0Z&version=1.1.0&request=GetMap&layers=fdi:GFDI_SFC&styles=&bbox=140.8879867553711,-39.76999832100885,151.05409088134766,-33.929998627233346&width=768&height=441&srs=EPSG:4326&format=application/openlayers#toggle

The project I am working on requires the NetCDF file to be updated on GeoServer on a daily basis.

When I was doing some testing with the 2017-12-28 data (3 days older), I was not able to preview the image for 2017-12-28T01:00:00.0Z

http://{GEOSERVER_DOMAIN}:8080/geoserver/fdi/wms?service=WMS&time=2017-12-28T01:00:00.0Z&version=1.1.0&request=GetMap&layers=fdi:GFDI_SFC&styles=&bbox=140.8879867553711,-39.76999832100885,151.05409088134766,-33.929998627233346&width=768&height=441&srs=EPSG:4326&format=application/openlayers#toggle

I was not able to preview data for any date and time earlier than 2018-01-02T01:00:00.0Z.

I have turned off Tile Caching and did disable and enable the NetCDF data store.

Is there anything else I need to look at/adjust to get this working? Is this a good practice to update the data source on a regular basis on GeoServer.

Edit - time dimension from the GetCapabilities request

<Dimension name="time" default="2017-12-27T18:00:00Z" units="ISO8601">
2017-12-27T18:00:00.000Z/2017-12-27T18:00:00.000Z/PT1S,2017-12-27T19:00:00.000Z/2017-12-27T19:00:00.000Z/PT1S,2017-12-27T20:00:00.000Z/2017-12-27T20:00:00.000Z/PT1S,2017-12-27T21:00:00.000Z/2017-12-27T21:00:00.000Z/PT1S,2017-12-27T22:00:00.000Z/2017-12-27T22:00:00.000Z/PT1S,2017-12-27T23:00:00.000Z/2017-12-27T23:00:00.000Z/PT1S,2017-12-28T00:00:00.000Z/2017-12-28T00:00:00.000Z/PT1S,2017-12-28T01:00:00.000Z/2017-12-28T01:00:00.000Z/PT1S,2017-12-28T02:00:00.000Z/2017-12-28T02:00:00.000Z/PT1S,2017-12-28T03:00:00.000Z/2017-12-28T03:00:00.000Z/PT1S,2017-12-28T04:00:00.000Z/2017-12-28T04:00:00.000Z/PT1S,
... ...
</Dimension>


texts: QGIS couldn't commit changes to layer on a database

texts:
i'm using QGIS on a database with a team,so all of our changes can be viewed at anytime
but sometimes we run into an issue sometimes when we save an error "couldn't commit changes" shows up and we have to exit the program and re-open it and start again from the last save point.. is there any particular reason why is this happening? any idea how to solve it ?


texts: How do I convert OS coordinates (X and Y) to longitude and latitude coordinates?

texts:
How do I convert OS coordinates (X and Y) - Eastings and Northings to longitude and latitude coordinates? For example X and Y below

X (Eastings): 347904

Y (Northings): 287484


texts: QGIS crashes when opening *.qgs file created with previous version?

texts:
I've upgraded QGIS to version 2.12.x (I tried both 2.12.2 and 2.12.4 versions) on Windows 7 x86-64.

QGIS crashes with the message 'minidump written to C:..." each time I try to load a .qgs file which has been created with a previous version.
However, there is no problem when I launch qgis-browser without opening the file.

The problem doesn't seem to come from the newer version, nor from the .qgis file because :

the .qgis file can be opened from an other computer running the version 2.12.2,
the .qgis file can be opened on the same computer but with an other Windows user account.

This led me to search the problem in my Windows user profile.
I tried to backup then remove all folders related to qgis in my user profile (.qgis, .qgis2, .matplotlib, .designer), but this has not solved the problem  :(

Is it possible that problem stems from my user registry ?

Assuming yes, what is the better way to solve it ?
Should I backup my NTUSER.DAT hive and then clean my user registry ?

Assuming yes, are there other keys than "HKEY_CURRENT_USERS\Software\QGIS" to delete ?